#########################################################
# R code to create ration and Upaij (aka predator foraging pref) for each size of each pred
# Kirstin Holsman
# Aug 2011
# Writes output to .dat file
########################################################
graphics.off()
rm(list=ls())					# clear R workspace
library(MASS)
library(lattice)
library(MCMCpack)
library(gam)
source("/Users/kkari/Documents/science/R_funKir/funKir.R") # load general functions
source("/Users/kkari/Documents/science/R_funKir/MSMfunctions.R") # load general functions
#########################################################
## 1. Read in the data
#########################################################
outfile<-"/Users/kkari/Documents/science/Projects/MSM/MSM12_modified/predwt/Wtpred2.dat"; app1=TRUE			# specific the .dat file to append (app1=TRUE) or not (app1=FALSE)
directorya<- "/Users/kkari/Documents/science/SVN/MSM/svn (trunk)/KIRM2"			# specify the working directory that contains the csv files
savefile<-"/Users/kkari/Documents/science/Projects/MSM/MSM12_modified/predwt"			# specify the folder to save the data files to
load(paste(directorya,"/msmdata2.RData",sep=""))
#LWdata2<-read.csv("/Users/kkari/Documents/science/Projects/MSM/msm_data/RACE_LenWt_out.csv")
Temp.data<-rbind(data.frame(Year=(1979:1981),EBS_bottomT=rep(mean(Temp.data$EBS_bottomT),3)),Temp.data,data.frame(Year=(2010:2011),EBS_bottomT=rep(mean(Temp.data$EBS_bottomT),2)))
overwrite<-1 # 0 means don't overwrite .dat file


LWAdata2<-read.csv("/Users/kkari/Documents/science/Projects/MSM/msm_data/LengthAtAge2.csv")
	dat<-LWAdata2[LWAdata2$GOAPOLL_PRED=="WALLEYE POLLOCK",]
	plk_WA<-na.omit(data.frame(age=dat$AGE,W=dat$WEIGHT,L=dat$LLENGTH/10))
	plk_WAT<-na.omit(data.frame(age=dat$AGE,W=dat$WEIGHT,L=dat$LLENGTH/10,T=dat$BOTTOMT))

	dat<-LWAdata2[LWAdata2$GOAPOLL_PRED=="PACIFIC COD",]
	pcod_WA<-na.omit(data.frame(age=dat$AGE,W=dat$WEIGHT,L=dat$LLENGTH/10))
	pcod_WAT<-na.omit(data.frame(age=dat$AGE,W=dat$WEIGHT,L=dat$LLENGTH/10,T=dat$BOTTOMT))

	dat<-LWAdata2[LWAdata2$GOAPOLL_PRED=="ARROWTOOTH FLOUNDR",]
	atf_WA<-na.omit(data.frame(age=dat$AGE,W=dat$WEIGHT,L=dat$LLENGTH/10))
	atf_WAT<-na.omit(data.frame(age=dat$AGE,W=dat$WEIGHT,L=dat$LLENGTH/10,T=dat$BOTTOMT))

LWdata2<-read.csv("/Users/kkari/Documents/science/Projects/MSM/msm_data/LWdata.csv")
	dat<-LWdata2[LWdata2$GOAPOLL_PRED=="WALLEYE POLLOCK",]
	LW_plk<-na.omit(data.frame(W=dat$WEIGHT,L=dat$LLENGTH/10))
	LW_plkT<-na.omit(data.frame(W=dat$WEIGHT,L=dat$LLENGTH/10,T=dat$BOTTOMT))
	dat<-LWdata2[LWdata2$GOAPOLL_PRED=="PACIFIC COD",]
	LW_pcod<-na.omit(data.frame(age=dat$AGE,W=dat$WEIGHT,L=dat$LLENGTH/10))
	LW_pcodT<-na.omit(data.frame(W=dat$WEIGHT,L=dat$LLENGTH/10,T=dat$BOTTOMT))

	dat<-LWdata2[LWdata2$GOAPOLL_PRED=="ARROWTOOTH FLOUNDR",]
	LW_atf<-na.omit(data.frame(age=dat$AGE,W=dat$WEIGHT,L=dat$LLENGTH/10))
	LW_atfT<-na.omit(data.frame(W=dat$WEIGHT,L=dat$LLENGTH/10,T=dat$BOTTOMT))

dat.use<-dat[dat$GOAPOLL_PRED=="WALLEYE POLLOCK",]

#setwd(directorya)															# set the working directory
#load("/Users/kkari/Documents/science/Projects/MSM/MSM_Rcode/msmdata.RData")
		#temps<-read.csv("temps.csv")							# temp data truncated to MSM years 1979:2002
	#	#Temp.data<-read.csv("BottomT.csv")						# read in temperature data
		#sub.lengths<-read.csv("/Users/kkari/Documents/science/Projects/MSM/msm_data/BS_AllLen_out.csv")					# read in the length data
		#sub.prey<-read.csv("/Users/kkari/Documents/science/Projects/MSM/msm_data/BS_allPrey_R.csv")						# read in the prey data
		#sub.prey<-data.frame(sub.prey)
		#prey<-rbind(sub.prey[which(sub.prey$Jmonth==6),],sub.prey[which(sub.prey$Jmonth==7),])
		#rm(sub.prey)
		#lengths<-rbind(sub.lengths[which(sub.lengths$MONTH==6),],sub.lengths[which(sub.lengths$MONTH==7),])
		#rm(sub.prey);rm(sub.lengths)
		#plk_WA<-read.csv("POLLOCK_SPECIMEN.csv")				# Pollock weight at age
		#pcod_WA<-read.csv("POLLOCK_SPECIMEN.csv")				# pcod weight at age - TOY DATA
		##arrow_WA<-read.csv("POLLOCK_SPECIMEN.csv")				# arrowtooth weight at age - TOY DATA
		##print("REMEMBER KIR TO LOAD PCOD AND ARROWTOOTH WEIGHT AT AGE DATA!!!")
		#LWpollock<-read.csv("LtoW_pollock.csv")
		#LWcod<-read.csv("LtoW_pcod.csv")
		#Temp.data2<-read.csv("BottomT.csv")						# read in temperature data
		#save(list=ls(),file="msmdata2.RData")
rm(LWpollock)
rm(LWcod)
prey$GearTemp[prey$GearTemp==-9]
#########################################################
## 1.2. Specify parameters and data below
#########################################################
model.parms<-rbind(
pollock=c(Ceq=2,Req=2,Weq=1,P=1,O2cal=13560,a=0.0143,CA=.380,CB=-.68,CQ=2.6,Tco=10,Tcm=15,RA=0.0075,RB=-.251,RQ=2.6,Tro=13,Trm=18,SDA=.125,Am=2,FA=.15,UA=.11),
pollock.juv=c(Ceq=2,Req=2,Weq=1,P=1,O2cal=13560,a=0.0143,CA=.340,CB=-.5875,CQ=6,Tco=8,Tcm=15,RA=0.0195,RB=-.26,RQ=4.6,Tro=15,Trm=18,SDA=.125,Am=1.4,FA=.2,UA=.11),
pollock.adult=c(Ceq=2,Req=2,Weq=1,P=1,O2cal=13560,a=0.0143,CA=.3,CB=-.5875,CQ=3.5,Tco=8,Tcm=15,RA=0.0137,RB=-.26,RQ=3.3,Tro=15,Trm=18,SDA=.125,Am=1.4,FA=.2,UA=.11),
pcod=c(Ceq=2,Req=2,Weq=1,P=1,O2cal=13560,a=0.0143,CA=0.03115127,CB=-.1342909,CQ=2.6,Tco=10,Tcm=15,RA=0.0075,RB=-.251,RQ=2.6,Tro=13,Trm=18,SDA=.125,Am=2,FA=.15,UA=.11),
arrowtooth=c(Ceq=1,Req=2,Weq=1,P=1,O2cal=13560,a=0.0143,CA=.0052,CB=-.22,CQ=.3763,Tco=15,Tcm=20,RA=0.0074,RB=-.496,RQ=3.077,Tro=13,Trm=18,SDA=.161,Am=1,FA=.15,UA=.1))
model.parms<-data.frame(model.parms)
n.spp<-3													# number of species, 1=pollock, 2=pcod, 3= arrowtooth
sp.dat<-list(spp=c("pollock","pcod","arrowtooth"),
	spp.p.pred=c(as.character(unique(prey$Species)[6]),as.character(unique(prey$Species)[3]),as.character(unique(prey$Species)[7])),		# pred names from prey file +_Wt, or _N, or _I
	spp.p.prey=as.character(c(names(prey)[45],names(prey)[129],"NA")),		# pred names from prey file
	spp.l=c("W. Pollock","P. Cod","Arrowtooth"),		# pred and prey names from the lengths file
	LW.a=rep(0,n.spp),										# LW a intercept for W=a*L^b		
	LW.b=rep(0,n.spp),
	CA=c(.380,0.03115127,.0052),#c(.380,.380,.0052),
	CB=c(-.68,-.1342909,-.22))	
names(sp.dat$LW.a)<-sp.dat$spp
names(sp.dat$LW.b)<-sp.dat$spp
sp.dat$LW.a2[3]<-5.682E-6 									# from Turnock GOA W = .003915 L^3.2232
sp.dat$LW.b2[3]<-3.1028 										# from Wilderbuer and Sample assessment (section 5) W=5.682E-6, b = 3.1028 (divide by 3)
for (i in 1:n.spp){
	eval(parse(text=paste(sp.dat$spp[i],"<-prey[prey$Species==sp.dat$spp.p.pred[i],]",sep="")))				#  diet wt pcod
	eval(parse(text=paste(sp.dat$spp[i],".l<-lengths[lengths$ECOPATH_PRED==sp.dat$spp.l[i],]",sep="")))		#  diet by lengths pcod.l
}
#########################################################
## 2. Calculate LW coeff for three species
#########################################################
par(mfrow=c(3,2))
dat<-LW_plk
dat2<-LW_plkT
LW.glm<-lm(log(W)~log(L),data=dat)
#LW.glm<-glm(W~L,data=dat, family=poisson(link = "log"))
LW.glmT<-lm(log(W)~log(L)*T,data=dat2)
sp.dat$LW.a[1]<-exp(coef(LW.glm)[1])
sp.dat$LW.b[1]<-coef(LW.glm)[2]
plot(dat$L,dat$W,pch=16,cex=.5,main="pollock")
points(dat2$L,exp(fitted(LW.glmT)),col=colors()[111],pch=16)
legend(200,3800,c("mean","mean+Temp"),col=c("red",colors()[111]),lty=1,lwd=2,box.lwd=0)

sp.dat$LW.a2[1]<-exp(coef(LW.glm)[1])
sp.dat$LW.b2[1]<-coef(LW.glm)[2]
LL<-1:120
predict.LW<-exp(predict.glm(LW.glm,newdata=data.frame(L=LL),type="response"))
#lines(LL,exp(coef(LW.glm)[1])*(LL^(coef(LW.glm)[2])),col="red", lwd=2)
lines(LL,predict.LW,col="red", lwd=2)
plot(density(dat$W),lwd=2)
lines(density(exp(fitted(LW.glm))),col="red")
lines(density(exp(fitted(LW.glmT))),col=colors()[111])
LW.glm.plk<-LW.glm
LW.glm.plkT<-LW.glmT

dat<-LW_pcod
dat2<-LW_pcodT
LW.glm<-lm(log(W)~log(L),data=dat)
LW.glmT<-lm(log(W)~log(L)*T,data=dat2)
sp.dat$LW.a[2]<-exp(coef(LW.glm)[1])
sp.dat$LW.b[2]<-coef(LW.glm)[2]
plot(dat$L,dat$W,pch=16,cex=.5,main="Pcod")
points(dat2$L,exp(fitted(LW.glmT)),col=colors()[111],pch=16)
sp.dat$LW.a2[2]<-exp(coef(LW.glm)[1])
sp.dat$LW.b2[2]<-coef(LW.glm)[2]
LL<-1:120
predict.LW<-exp(predict.glm(LW.glm,newdata=data.frame(L=LL),type="response"))
lines(LL,exp(coef(LW.glm)[1])*(LL^(coef(LW.glm)[2])),col="red", lwd=2)
lines(LL,predict.LW,col="red", lwd=2)
plot(density(dat$W),lwd=2)
lines(density(exp(coef(LW.glm)[1])*(dat$L^(coef(LW.glm)[2]))),col="red")
lines(density(exp(fitted(LW.glmT))),col=colors()[111])
LW.glm.pcod<-LW.glm
LW.glm.pcodT<-LW.glmT

dat<-LW_atf
dat2<-LW_atfT
LW.glm<-lm(log(W)~log(L),data=dat)
LW.glmT<-lm(log(W)~log(L)*T,data=dat2)
sp.dat$LW.a[3]<-exp(coef(LW.glm)[1])
sp.dat$LW.b[3]<-coef(LW.glm)[2]
plot(dat$L,dat$W,pch=16,cex=.5,main="Arrowtooth")
points(dat2$L,exp(fitted(LW.glmT)),col=colors()[111],pch=16)
#sp.dat$LW.a2[3]<-exp(coef(LW.glm)[1])
#sp.dat$LW.b2[3]<-coef(LW.glm)[2]
LL<-1:120
predict.LW<-exp(predict.glm(LW.glm,newdata=data.frame(L=LL),type="response"))
lines(LL,exp(coef(LW.glm)[1])*(LL^(coef(LW.glm)[2])),col="red", lwd=2)
lines(LL,predict.LW,col="red", lwd=2)
lines(LL,5.682E-6*((LL*10)^3.1028),col="green")
plot(density(dat$W))
lines(density(exp(coef(LW.glm)[1])*(dat$L^(coef(LW.glm)[2]))),col="red")
lines(density(exp(fitted(LW.glmT))),col=colors()[111])
LW.glm.atf<-LW.glm
LW.glm.atfT<-LW.glmT

prey$plkWt.ratio<-prey$WALLEYE.POLLOCK_Wt/prey$TotWt
prey$pcodWt.ratio<-prey$PACIFIC.COD_Wt/prey$TotWt
prey$otherWt.ratio<-(prey$TotWt-prey$PACIFIC.COD_Wt-prey$WALLEYE.POLLOCK_Wt)/prey$TotWt

#########################################################
## 2. Calculate Upaij for each pred - eg. test.21$Upaij is the Upaij for pcod eating pollock
#########################################################
l.bins<-seq(0,140,1)
ppmat.p<-matrix(0,length(l.bins),length(l.bins))
colnames(ppmat.p)<-l.bins
rownames(ppmat.p)<-l.bins
par(mfrow=c(3,1))
for (predd in 1:n.spp){	
	eval(parse(text=paste("sub.dat<-",sp.dat$spp[predd],".l",sep="")))
	eval(parse(text=paste("sub.dat2<-",sp.dat$spp[predd],sep="")))
	sub.dat2$Wt<-sp.dat$LW.a[predd]*(sub.dat2$Length^sp.dat$LW.b[predd])	
	eval(parse(text=paste(sp.dat$spp[predd],"$Wt<-sub.dat2$Wt",sep="")))
	eval(parse(text=paste("S.",predd,"<-data.frame(S=tapply(sub.dat2$TotWt/sub.dat2$Wt,sub.dat2$Length,mean.na))",sep=""))) # mean stomach weights
	eval(parse(text=paste("S.",predd,"$S.sd<-tapply(sub.dat2$TotWt/sub.dat2$Wt,sub.dat2$Length,sd.na)",sep="")))
	eval(parse(text=paste("S.",predd,"$S.se<-tapply(sub.dat2$TotWt/sub.dat2$Wt,sub.dat2$Length,se.na)",sep="")))
	eval(parse(text=paste("S.",predd,"$C<-(tapply((sub.dat2$TotWt/sub.dat2$Wt)*(24*0.0134*exp(0.115*sub.dat2$GearTemp)),sub.dat2$Length,mean.na))",sep="")))
	eval(parse(text=paste("S.",predd,"$C.sd<-(tapply((sub.dat2$TotWt/sub.dat2$Wt)*(24*0.0134*exp(0.115*sub.dat2$GearTemp)),sub.dat2$Length,sd.na))",sep="")))
	eval(parse(text=paste("S.",predd,"$C.se<-(tapply((sub.dat2$TotWt/sub.dat2$Wt)*(24*0.0134*exp(0.115*sub.dat2$GearTemp)),sub.dat2$Length,se.na))",sep="")))
	eval(parse(text=paste("S.",predd,"$Length<-as.numeric(rownames(S.",predd,"))",sep="")))
	eval(parse(text=paste("S.",predd,"$Wt<-sp.dat$LW.a[predd]*(S.",predd,"$Length^sp.dat$LW.b[predd])",sep="")))
	np<-c(1,4,5)
	FT<-bioenergetics(par=model.parms[np[predd],],W=S.1$Wt[1],TempC=sub.dat2$GearTemp,P=1,Eprey=3000,Efish=3000)$fTc
	eval(parse(text=paste(sp.dat$spp[predd],"$fT<-FT",sep="")))
	eval(parse(text=paste("sub.dat2<-",sp.dat$spp[predd],sep="")))
	eval(parse(text=paste("S.",predd,"$C.hat<-(tapply((sp.dat$CA[predd]*(sub.dat2$Wt^sp.dat$CB[predd])*sub.dat2$fT),sub.dat2$Length,mean.na))",sep="")))
	eval(parse(text=paste("S.",predd,"$C.hatTot<-(tapply((sp.dat$CA[predd]*(sub.dat2$Wt^sp.dat$CB[predd])*sub.dat2$fT*sub.dat2$Wt),sub.dat2$Length,mean.na))",sep="")))
	
	eval(parse(text=paste("S.",predd,"$C.hat.noFT<-sp.dat$CA[predd]*(S.",predd,"$Wt^sp.dat$CB[predd])",sep="")))
	eval(parse(text=paste("ST.",predd,"<-data.frame(S=tapply(sub.dat2$TotWt/sub.dat2$Wt,sub.dat2$GearTemp,max.na))",sep="")))
	eval(parse(text=paste("ST.",predd,"$T<-as.numeric(rownames(ST.",predd,"))",sep="")))
	eval(parse(text=paste("S.",predd,"$fT<-(tapply(sub.dat2$fT,sub.dat2$Length,mean.na))",sep="")))
	eval(parse(text=paste("S.",predd,"$meanT<-(tapply(sub.dat2$GearTemp,sub.dat2$Length,mean.na))",sep="")))

#	model.b<-lm(C~L+L2,data=list(C=sub.dat2$TotWt/(sub.dat2$Wt),L=sub.dat2$Length,L2=sub.dat2$Length^2))
	tt<-tapply(sub.dat2$TotWt/(sub.dat2$Wt),sub.dat2$Length,mean.na)
	plot(as.numeric(names(tt)),tt,pch=16)
	model.b<-lm(C~L+I(L^2)+I(L^3)+I(L^4),data=list(C=sub.dat2$TotWt/(sub.dat2$Wt),L=sub.dat2$Length))
	
	tt2<-predict(model.b,newdata=list(L=as.numeric(names(tt))))

		lines(as.numeric(names(tt)),tt2,pch=16,col="red",lwd=2)
	#model.b<-lm(logC~logW,data=list(logC=log(sub.dat2$TotWt/sub.dat2$Wt),logW=log(sub.dat2$Wt),L2=sub.dat2$Length^2))
	#model.b2<-glm(C~L,data=list(C=0.00001+(sub.dat2$TotWt/sub.dat2$Wt),L=sub.dat2$Length,L2=sub.dat2$Length^2),family=Gamma(link=log))
	sp.dat$aL[predd]<-(coef(model.b)[1])
	sp.dat$bL[predd]<-coef(model.b)[2]
	sp.dat$bL2[predd]<-coef(model.b)[3]
	sp.dat$bL3[predd]<-coef(model.b)[4]
	sp.dat$bL4[predd]<-coef(model.b)[5]
	sp.dat$bL5[predd]<-coef(model.b)[6]
	
	eval(parse(text=paste("S.",predd,"$S.MSM<-(tapply((fitted(model.b)),sub.dat2$Length,mean.na))",sep="")))
	eval(parse(text=paste("S.",predd,"$C.MSM<-(tapply((fitted(model.b))*(24*0.0134*exp(0.115*sub.dat2$GearTemp)),sub.dat2$Length,mean.na))",sep="")))
	eval(parse(text=paste("S.",predd,"$C.MSM.ft<-(tapply((fitted(model.b))*FT,sub.dat2$Length,mean.na))",sep="")))
	for (preyy in 1:(n.spp-1)){
		sub.dat<-sub.dat[sub.dat$ECOPATH_PREY==sp.dat$spp.l[preyy],]
		sp.t<-strsplit(sp.dat$spp.p.prey[preyy],"_")[[1]][1]
		eval(parse(text=paste("test.",predd,preyy,"<-stom.wt(sub.dat=sub.dat,a.prey=sp.dat$LW.a[preyy],b.prey=sp.dat$LW.b[preyy],pred.spp=sp.dat$spp[predd],prey.spp=sp.t)",sep="")))
	}
	rm(sub.dat);rm(sub.dat2)
}
ST.1<-ST.1[ST.1$S<=.2,]
ST.2<-ST.2[ST.2$S<=.5,]
S.1$S.hat2<-model.parms$CA[3]*(S.1$Wt^model.parms$CB[3])
S.1$S.hat2[S.1$Length<20]<-model.parms$CA[2]*(S.1$Wt[S.1$Length<20]^model.parms$CB[2])

rm(test)

stomWT=matrix(0,3,length(1:150));stomWT=data.frame(stomWT)
colnames(stomWT)<-1:150
LL<-as.numeric(colnames(stomWT))
sub.dat1<-pollock
sub.dat2<-pollock[pollock$TotWt>0.001,]
sub.dat2<-sub.dat2[sub.dat2$Length<80,]
sub.dat2<-sub.dat2[sub.dat2$Length>5,]
meanS<-tapply(sub.dat2$TotWt/sub.dat2$Wt,sub.dat2$Length,mean.na)
meanS.aka<-tapply(sub.dat1$TotWt/sub.dat1$Wt,sub.dat1$Length,mean.na)
plot(as.numeric(names(meanS)),meanS,ylim=c(0,.1))
points(as.numeric(names(meanS.aka)),meanS.aka,ylim=c(0,.1))

model.c<-lm(logC~logW,data=list(logC=log(sub.dat2$TotWt/sub.dat2$Wt),logW=log(sub.dat2$Wt),L2=sub.dat2$Length^2))
points(sub.dat2$Length,exp(fitted(model.c)),pch=16)
aa<-exp(coef(model.c)[1])
bb<-coef(model.c)[2]
WW<-sp.dat$LW.a[1]*LL^sp.dat$LW.a[2]
stomWT[1,]<-aa*WW^bb
lines(LL,stomWT[1,])



plot(LL,stomWT[1,]);lines(S.1$Length,S.1$S)
############################################################
## Plot Temperature as mean temp, and biomass or numerically weighted temperature for each spp
############################################################
colls<-c(336,88,148)
dev.new(width=6,height=6.5)
	layout(rbind(c(1,0),
				c(2,0))
		,widths=c(1,0.01),
		heights=rep(1,2))
	par(mar=c(4,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(1,1,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(2)

#Temp biomass by year
temp.add<-tapply(pollock$Wt*pollock$GearTemp,pollock$Year,sum.na)/tapply(pollock$Wt,pollock$Year,sum.na)
Temp.data$pollock.T<-c(rep(mean.na(temp.add),3),temp.add,rep(mean.na(temp.add),2))
tt<-tapply(pcod$Wt*pcod$GearTemp,pcod$Year,sum.na)/tapply(pcod$Wt,pcod$Year,sum.na)
Temp.data$pcod.T<- rep(mean(tt),length(Temp.data$Year))
Temp.data$pcod.T[na.omit(match(as.numeric(names(tt)),Temp.data$Year))]<-tt[na.omit(match(Temp.data$Year,as.numeric(names(tt))))]
tt<-tapply(arrowtooth$Wt*arrowtooth$GearTemp,arrowtooth$Year,sum.na)/tapply(arrowtooth$Wt,arrowtooth$Year,sum.na)
Temp.data$arrowtooth.T<- rep(mean(tt),length(Temp.data$Year))
Temp.data$arrowtooth.T[na.omit(match(as.numeric(names(tt)),Temp.data$Year))]<-tt[na.omit(match(Temp.data$Year,as.numeric(names(tt))))]

plot(Temp.data$Year,Temp.data[,2],type="l",lwd=2,lty=2,ylim=c(0,4),xlab="",main="Biomass weighted Temp",axes=FALSE);axis(1);axis(2)
points(Temp.data$Year,Temp.data$pollock.T,type="l",col="red",lwd=2)
points(Temp.data$Year,Temp.data$pcod.T,type="l",col="blue",lwd=2)
points(Temp.data$Year,Temp.data$arrowtooth.T,type="l",col="green4",lwd=2)
legend(2005,4,c("Temp index","Pollock","PCod","Arrowtooth"),col=c("black","red","blue","green4"),lty=c(2,1,1,1),lwd=2,box.col=FALSE,cex=.7)
# Temp numerically

temp.add<-tapply(pollock$GearTemp,pollock$Year,mean.na)
Temp.data$pollock.T.num<-c(rep(mean.na(temp.add),3),temp.add,rep(mean.na(temp.add),2))
tt<-tapply(pcod$GearTemp,pcod$Year,mean.na)
Temp.data$pcod.T.num<- rep(mean(tt),length(Temp.data$Year))
Temp.data$pcod.T.num[na.omit(match(as.numeric(names(tt)),Temp.data$Year))]<-tt[na.omit(match(Temp.data$Year,as.numeric(names(tt))))]
tt<-tapply(arrowtooth$GearTemp,arrowtooth$Year,mean.na)
Temp.data$arrowtooth.T.num<- rep(mean(tt),length(Temp.data$Year))
Temp.data$arrowtooth.T.num[na.omit(match(as.numeric(names(tt)),Temp.data$Year))]<-tt[na.omit(match(Temp.data$Year,as.numeric(names(tt))))]

plot(Temp.data$Year,Temp.data[,2],type="l",lwd=2,lty=2,ylim=c(0,4),xlab="Year",main="Numerically weighted Temp",axes=FALSE);axis(1);axis(2)
points(Temp.data$Year,Temp.data$pollock.T.num,type="l",col="red",lwd=2)
points(Temp.data$Year,Temp.data$pcod.T.num,type="l",col="blue",lwd=2)
points(Temp.data$Year,Temp.data$arrowtooth.T.num,type="l",col="green4",lwd=2)



#########################################################
## Dirichlet MLE 
#########################################################
colls<-c(336,88,148)
dev.new(width=6,height=6.5)
	layout(rbind(c(1,0,2),
				c(3,0,4),
				c(5,0,6))
		,widths=c(1,0.01,1),
		heights=rep(1,3))
	par(mar=c(4,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(1,1,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(6)

##### plot diet comp as a function of length of the predator
dat2<-pollock
#dat1.y<-data.frame(Length=dat2$Length,Wt=dat2$Wt,plk=dat2$WALLEYE.POLLOCK_Wt,pcod=dat2$PACIFIC.COD_Wt,Total=dat2$TotWt)
dat1.y<-data.frame(Length=dat2$Length,Wt=dat2$Wt,plk=dat2$WALLEYE.POLLOCK_Wt/dat2$Wt,pcod=dat2$PACIFIC.COD_Wt/dat2$Wt,Total=dat2$TotWt/dat2$Wt)
dat1.y$other<-dat1.y$Total-(dat1.y$plk+dat1.y$pcod)
y<-tapply(dat1.y$plk,dat1.y$Length,mean.na);x<-as.numeric(names(y))
plot(x,y,type="l",axes=FALSE,xlim=c(1,max(x)),col=colors()[colls[2]+1]);axis(1);axis(2)
y<-tapply(dat1.y$pcod,dat1.y$Length,mean.na);x<-as.numeric(names(y))
lines(x,y,type="l",col=colors()[colls[3]+1])
y<-tapply(dat1.y$plk/dat1.y$Total,dat1.y$Length,mean.na);x<-as.numeric(names(y))
nobs<-length(y)
mtext(bquote(bold("Wt of prey (as prop of BW)")),cex=.7,side=3,line=-2,outer=FALSE)
plot(x,y,type="l",axes=FALSE,ylim=c(0,1),xlim=c(1,max(x)));axis(1);axis(2)
polygon(c(0,max(x),max(x),0),c(1,1,0,0),col=colors()[colls[1]],border=FALSE)
y<-tapply((dat1.y$pcod+dat1.y$plk)/dat1.y$Total,dat1.y$Length,mean.na)
polygon(c(x,x[nobs:1]),c(y,rep(0,nobs)),col=colors()[colls[2]],border=colors()[colls[2]+1])
y<-tapply((dat1.y$pcod)/dat1.y$Total,dat1.y$Length,mean.na)
polygon(c(x,x[nobs:1]),c(y,rep(0,nobs)),col=colors()[colls[3]],border=colors()[colls[3]+2])
mtext(bquote(bold("Pred: Pollock")),cex=.9,side=3,line=-1,outer=TRUE)
mtext(bquote(bold("Proportion of Stom")),cex=.7,side=3,line=-2,outer=FALSE)

##### plot diet comp as a function of length of the predator
dat2<-pcod
#dat1.y<-data.frame(Length=dat2$Length,Wt=dat2$Wt,plk=dat2$WALLEYE.POLLOCK_Wt,pcod=dat2$PACIFIC.COD_Wt,Total=dat2$TotWt)
dat1.y<-data.frame(Length=dat2$Length,Wt=dat2$Wt,plk=dat2$WALLEYE.POLLOCK_Wt/dat2$Wt,pcod=dat2$PACIFIC.COD_Wt/dat2$Wt,Total=dat2$TotWt/dat2$Wt)
dat1.y$other<-dat1.y$Total-(dat1.y$plk+dat1.y$pcod)
y<-tapply(dat1.y$plk,dat1.y$Length,mean.na);x<-as.numeric(names(y))
plot(x,y,type="l",axes=FALSE,xlim=c(1,max(x)),col=colors()[colls[2]+1]);axis(1);axis(2)
y<-tapply(dat1.y$pcod,dat1.y$Length,mean.na);x<-as.numeric(names(y))
lines(x,y,type="l",col=colors()[colls[3]+1])
y<-tapply(dat1.y$plk/dat1.y$Total,dat1.y$Length,mean.na);x<-as.numeric(names(y))
nobs<-length(y)
#mtext(bquote(bold("Wt of prey (as prop of BW)")),cex=.7,side=3,line=-2,outer=FALSE)
plot(x,y,type="l",axes=FALSE,ylim=c(0,1),xlim=c(1,max(x)));axis(1);axis(2)
polygon(c(0,max(x),max(x),0),c(1,1,0,0),col=colors()[colls[1]],border=FALSE)
y<-tapply((dat1.y$pcod+dat1.y$plk)/dat1.y$Total,dat1.y$Length,mean.na)
polygon(c(x,x[nobs:1]),c(y,rep(0,nobs)),col=colors()[colls[2]],border=colors()[colls[2]+1])
y<-tapply((dat1.y$pcod)/dat1.y$Total,dat1.y$Length,mean.na)
polygon(c(x,x[nobs:1]),c(y,rep(0,nobs)),col=colors()[colls[3]],border=colors()[colls[3]+2])
mtext(bquote(bold("Pred: Pcod")),cex=.9,side=3,line=-15,outer=TRUE)
#mtext(bquote(bold("Proportion of Stom")),cex=.7,side=3,line=-2,outer=FALSE)
text(20,.8,bquote(bold("Other prey")),col=colors()[colls[3]+10])
text(90,.2,bquote(bold("Pollock")),col=colors()[497])
text(125,.7,bquote(bold("Pcod")),col=colors()[75])

##### plot diet comp as a function of length of the predator
dat2<-arrowtooth
#dat1.y<-data.frame(Length=dat2$Length,Wt=dat2$Wt,plk=dat2$WALLEYE.POLLOCK_Wt,pcod=dat2$PACIFIC.COD_Wt,Total=dat2$TotWt)
dat1.y<-data.frame(Length=dat2$Length,Wt=dat2$Wt,plk=dat2$WALLEYE.POLLOCK_Wt/dat2$Wt,pcod=dat2$PACIFIC.COD_Wt/dat2$Wt,Total=dat2$TotWt/dat2$Wt)
dat1.y$other<-dat1.y$Total-(dat1.y$plk+dat1.y$pcod)
y<-tapply(dat1.y$plk,dat1.y$Length,mean.na);x<-as.numeric(names(y))
plot(x,y,type="l",axes=FALSE,xlim=c(1,max(x)),col=colors()[colls[2]+1]);axis(1);axis(2)
y<-tapply(dat1.y$pcod,dat1.y$Length,mean.na);x<-as.numeric(names(y))
lines(x,y,type="l",col=colors()[colls[3]+1])
y<-tapply(dat1.y$plk/dat1.y$Total,dat1.y$Length,mean.na);x<-as.numeric(names(y))
nobs<-length(y)
#mtext(bquote(bold("Wt of prey (as prop of BW)")),cex=.7,side=3,line=-2,outer=FALSE)
plot(x,y,type="l",axes=FALSE,ylim=c(0,1),xlim=c(1,max(x)));axis(1);axis(2)
polygon(c(0,max(x),max(x),0),c(1,1,0,0),col=colors()[colls[1]],border=FALSE)
y<-tapply((dat1.y$pcod+dat1.y$plk)/dat1.y$Total,dat1.y$Length,mean.na)
polygon(c(x,x[nobs:1]),c(y,rep(0,nobs)),col=colors()[colls[2]],border=colors()[colls[2]+1])
y<-tapply((dat1.y$pcod)/dat1.y$Total,dat1.y$Length,mean.na)
polygon(c(x,x[nobs:1]),c(y,rep(0,nobs)),col=colors()[colls[3]],border=colors()[colls[3]+2])
mtext(bquote(bold("Pred: Arrowtooth")),cex=.9,side=3,line=-32,outer=TRUE)
#mtext(bquote(bold("Proportion of Stom")),cex=.7,side=3,line=-2,outer=FALSE)
dev.print(device=postscript, paste("K.eps",sep=""), onefile=FALSE, horizontal=FALSE)




skip<-1
if(skip==0){
	# simulate data and test with mle
	LL<-seq(1,100,1)
	nobs<-length(LL)
	a<-c(-80,0.5,20)
	b<-c(1.3,1.1,-2.4)
	m<-100
	alpha.sim<-matrix(0,nobs,3)
	Y.sim2<-matrix(0,nobs,3)
	for (i in 1:nobs){
		alpha.sim[i,]<-inv.logit(a+LL[i]*b)*m
		#alpha.sim[i,]<-(a+(LL[i]*b))		# larger the alpha value the tighter the disribution
		Y.sim2[i,]<-(rdirichlet(1,alpha.sim[i,]))
		if(any(Y.sim2[i,]==0)){
				Y.sim2[i,which(Y.sim2[i,]==0)]<-10e-20
				Y.sim2[i,]<-Y.sim2[i,]/sum(Y.sim2[i,])
		}
	}
	dat2.sim<-data.frame(Length=LL,Y.sim2)
	par(mfrow=c(3,1))
	plot(dat2.sim[,2],ylim=c(0,1),pch=16)
	plot(dat2.sim[,3],ylim=c(0,1),pch=16)
	plot(dat2.sim[,4],ylim=c(0,1),pch=16)
	## function for dirichlet MLE
	dir.mle = function(par,data) {
		data<-data
		par<-par
		datt<-data$dat1
		k<-data$k				# number of prey categories
		Y<-datt[,2:4]
		nobs<-length(Y[,1])
		Length<-datt[,1]
		a<-(par[1:3]) 		
		b<-(par[4:6])		# b must be >0
		m<-exp(par[7])
		ll<-rep(0,nobs)
		alpha.hat<-matrix(0,nobs,3)
		for(i in 1:nobs){
			alpha.hat[i,]<-inv.logit(a+(Length[i])*b)*m
			ll[i]<-ddirichlet(x=Y[i,],alpha=alpha.hat[i,])
			## remove 0's since dirichlet can't handle them - ok to do?
			if(any(Y[i,]==0)){
				Y[i,which(Y[i,]==0)]<-10e-20
				Y[i,]<-Y[i,]/sum(Y[i,])
			}		
			ll[i]<-ddirichlet(x=Y[i,],alpha=alpha.hat[i,])
			ll[i]<-ifelse(ll[i]==0,1e-200,ll[i])  # penalize the LL for values of 0 (otherwise log(0)<--Inf)
		}
		-sum(log(ll))
	}
	dir.mle(par=c(a=(c(-11,-1,1)),logb=log(c(1,1,1)),logm=log(10)),data=list(dat1=dat2.sim,k=3) )			# test the function to be sure it returns NLL
	mm2<-optim(fn=dir.mle,par=c(a=(c(-10,-1,1)),logb=log(c(1,1,1)),logm=log(100)),data=list(dat1=dat2.sim,k=3),control=list(trace=1,maxit=1000),hessian=TRUE)  # now use mle to find parms
	a.hat2<-(mm2$par[1:3])
	b.hat2<-exp(mm2$par[4:6])
	m.hat2<-exp(mm2$par[7])
	round(cbind(a,a.hat2,b,b.hat2,m,m.hat2),4)
	## Now calculate Y.hat from a.hat, b.hat
	alpha.hat<-matrix(0,nobs,3)			# pre-allocate
	Y.hat2<-alpha.hat					# pre-allocate
	for (i in 1:nobs){
		alpha.hat[i,]<-inv.logit(a.hat2+(dat2.sim[i,1]*b.hat2))*m.hat2
		Y.hat2[i,]<-alpha.hat[i,]/sum(alpha.hat[i,])
	}	
	plot(dat2.sim[,1],dat2.sim[,2],ylim=c(0,1),pch=16)
	lines(Y.hat2[,1],col="red",lwd=2)
	plot(dat2.sim[,1],dat2.sim[,3],ylim=c(0,1),pch=16)
	lines(Y.hat2[,2],col="red",lwd=2)
	plot(dat2.sim[,1],dat2.sim[,4],ylim=c(0,1),pch=16)
	lines(Y.hat2[,3],col="red",lwd=2)
	plot(alpha.hat[,1],ylim=c(0,1))
	plot(alpha.hat[,2],ylim=c(0,1))
	plot(alpha.hat[,3],ylim=c(0,1))
	
	#######################################
	### plot K ~ Length for each pred
	#######################################	
	dev.new(width=6,height=6.5)
		layout(rbind(c(1,0,4,0,7),
					c(2,0,5,0,8),
					c(3,0,6,0,9))
			,widths=c(1,0.01,1,0.01,1),
			heights=rep(1,3))
		par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
		par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
		par(oma=c(5,5,1,0))# outer margins of graph (bottom,left, top, right)
		layout.show(9)	
	sp.dat$a.hat<-matrix(0,n.spp,n.spp)
	sp.dat$b.hat<-matrix(0,n.spp,n.spp)	
	for (predd in 1:3){
		eval(parse(text=paste("dat2<-",sp.dat$spp[predd],sep="")))
		dat1<-data.frame(Length=dat2$Length,plk.prop=dat2$WALLEYE.POLLOCK_Wt/dat2$TotWt,pcod.prop=dat2$PACIFIC.COD_Wt/dat2$TotWt)
		dat1$other.prop<-1-(dat1$plk.prop+dat1$pcod.prop)
		dat1<-dat1[-which(is.na(rowSums(dat1[,2:4]))),]
		# subset of data
		dat1.use<-data.frame(Length=as.numeric(names(tapply(dat1[,2],dat1[,1],mean))),plk.prop=tapply(dat1[,2],dat1[,1],mean),pcod.prop=tapply(dat1[,3],dat1[,1],mean))
		dat1.use$other.prop<-1-(dat1.use$plk.prop+dat1.use$pcod.prop)
		dir.mle(par=c(loga=log(c(.2,.01,2)),logb=log(c(2,1,.4))),data=list(dat1=dat1.use,k=3))
		mm<-optim(fn=dir.mle,par=c(loga=log(c(.2,.01,2)),logb=log(c(2,1,.4))),data=list(dat1=dat1.use,k=3),control=list(trace=0,maxit=5000))
		nobs<-length(dat1.use[,1])
		a.hat<-exp(mm$par[1:3])
		b.hat<-exp(mm$par[4:6])
		sp.dat$a.hat[predd,]<-a.hat
		sp.dat$b.hat[predd,]<-b.hat		
		alpha.hat<-matrix(0,nobs,3)			# pre-allocate
		Y.hat<-alpha.hat					# pre-allocate
		
		for (i in 1:nobs){
			alpha.hat[i,]<-(a.hat+((LL[i])*b.hat))
			Y.hat[i,]<-alpha.hat[i,]/sum(alpha.hat[i,])
		}
		plot(dat1.use[,2],ylim=c(0,1),pch=16,main=sp.dat$spp[predd])
		lines(Y.hat[,1],col="red")
		plot(dat1.use[,3],ylim=c(0,1),pch=16)
		lines(Y.hat[,2],col="red")
		plot(dat1.use[,4],ylim=c(0,1),pch=16)
		lines(Y.hat[,3],col="red")		
	}
	dev.print(device=postscript, paste("dirFits.eps",sep=""), onefile=FALSE, horizontal=FALSE)

}#end skip	

#########################################################
## Gamma dist for K
#########################################################




sp.dat$shape<-matrix(0,n.spp,n.spp)
rownames(sp.dat$shape)<-sp.dat$spp
colnames(sp.dat$shape)<-sp.dat$spp
sp.dat$scale<-sp.dat$shape
sp.dat$a.glm<-sp.dat$shape
sp.dat$b.glm<-sp.dat$shape
dev.new(width=6.5,height=5)
	layout(rbind(c(1,0,3,0,5),
				c(2,0,4,0,6))
		,widths=c(1,0.01,1,0.01,1),
		heights=rep(1,2))
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(5,5,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(6)
## derive parameters for the gamma function for proportion of each prey item in each pred diet
maxK<-matrix(1,n.spp,n.spp)
maxKL<-maxK
predd<-1
preyy<-1

for (predd in 1:n.spp){
	for (preyy in 1:(n.spp-1)){
		eval(parse(text=paste("dat.use<-test.",predd,preyy,"$wtot.ratio.mean",sep="")))
		dat1<-na.omit(data.frame(Length=as.numeric(names(dat.use)),prop=dat.use))
		dat1<-na.omit(data.framet(Length=prey.
		dat1$prop[dat1$prop==0]<-0.0001
		if(preyy>1){
			dat1<-dat1[dat1$prop<0.04,]			
		}
		m3<-glm(prop~length,data=list(length=(dat1$Length),prop=dat1$prop),family=Gamma(link="log"))
		maxK[predd,preyy]<-max(dat1$prop)
		tl<-dat1$Length[match(maxK[predd,preyy],dat1$prop)]
		maxKL[predd,preyy]<-tl[1]
		shape<-as.numeric(gamma.shape(m3)[1])
		sp.dat$shape[predd,preyy]<-as.numeric(gamma.shape(m3)[1])
		scalee<-as.numeric(mean(dat1$prop)/shape)
		sp.dat$scale[predd,preyy]<-as.numeric(mean(dat1$prop)/shape)
		sp.dat$a.glm[predd,preyy]<-(coef(m3)[1])
		sp.dat$b.glm[predd,preyy]<-(coef(m3)[2])	
		plot(dat1$Length,dat1$prop,axes=FALSE,main=paste(sp.dat$spp[predd],"-->",sp.dat$spp[preyy]),ylim=c(0,max(dat1$prop)*1.1),pch=16);axis(1);axis(2)
		#points(dat1$Length,fitted(m3),pch=16,col="red")
		lines(dat1$Length,exp(coef(m3)[1]+dat1$Length*coef(m3)[2]),col="red",lwd=2)
		Kprop<-exp(coef(m3)[1]+dat1$Length*coef(m3)[2])
		Kprop[Kprop>maxK[predd,preyy]]<-maxK[predd,preyy]
		if(preyy==1){Kprop[Kprop>maxK[predd,preyy]]<-0}
		dat2<-dat1[dat1$Length<=maxKL[predd,preyy],]
		dat2<-dat1
		
		model.b<-glm(P~L+I(L^2)+I(L^3)+I(L^4)+I(L^5)+I(L^6),data=data.frame(P=dat2$prop,L=dat2$Length),family=Gamma(link="log"))

		gam.1<-gam(dat2$prop~s(dat2$Length,1),family=Gamma(link="log"))
		gam.2<-gam(dat2$prop~s(dat2$Length,2),family=Gamma(link="log"))
		gam.3<-gam(dat2$prop~s(dat2$Length,3),family=Gamma(link="log"))
		gam.4<-gam(dat2$prop~s(dat2$Length,4),family=Gamma(link="log"))
		gam.5<-gam(dat2$prop~s(dat2$Length,5),family=Gamma(link="log"))
		gam.6<-gam(dat2$prop~s(dat2$Length,6),family=Gamma(link="log"))
			gam.7<-gam(dat2$prop~s(dat2$Length,7),family=Gamma(link="log"))
				gam.8<-gam(dat2$prop~s(dat2$Length,8),family=Gamma(link="log"))
		aa<-AIC(gam.1,gam.2,gam.3,gam.4,gam.5,gam.6,gam.7,gam.8)
		aa$num<-1:length(aa)
		aa$deltaAIC<-aa$AIC-min(aa$AIC)
		aa$rank<-rank(aa$deltaAIC)
		aa<-aa[order(aa$rank),]
		eval(parse(text=paste("gam.use<-",rownames(aa)[1],sep="")))
		tt2<-(predict(gam.use,se.fit = TRUE))
		tt2.m<-exp(tt2$fit)
		tt2.u<-exp(tt2$fit+1.95*tt2$se.fit)
		tt2.l<-exp(tt2$fit-1.95*tt2$se.fit)
		mtext(side=3,rownames(aa)[1],line=-1)
		#lines(dat1$Length,Kprop,col="blue",lwd=2)
		lines(dat2$Length,tt2.m,col="black",lwd=2)
		lines(dat2$Length,tt2.u,col="green2",lwd=2,lty=2)
		lines(dat2$Length,tt2.l,col="green2",lwd=2,lty=2)
		
		points(dat1$Length,dat1$prop,pch=16)
		m3<-glm(prop~length,data=list(length=(dat1$Length),prop=dat1$prop),family=Gamma(link="log"))
		#points(dat1$Length,fitted(m3),pch=16,col="red")
		#lines(dat1$Length,exp(coef(m3)[1])*exp(dat1$Length*coef(m3)[2]))
		rate<-1/scalee
		disp<-1/shape
	
	}
}
sp.dat$maxK<-maxK
# write smoothed K to .dat file
#########################################################
## NOW predict Upaij.hat from st.biomass*switch*Kpij
#########################################################
sp.dat$lim.glm.a<-rep(0,n.spp)
#rownames(sp.dat$lim.glm)<-sp.dat$spp
names(sp.dat$lim.glm.a)<-sp.dat$spp
sp.dat$lim.glm.b<-sp.dat$lim.glm.a

par(mfrow=c(3,2))
for (predd in 1:3){
	eval(parse(text=paste("subb<-",sp.dat$spp[predd],sep="")))
	eval(parse(text=paste("subb.l<-",sp.dat$spp[predd],".l",sep="")))
	max.prey.size<-tapply(subb.l$PREY_SIZE_CM,subb.l$PRED_LEN,max)		
	max.prey.size<-max.prey.size[max.prey.size>0]
	max.prey.size<-max.prey.size[as.numeric(names(max.prey.size))<80]
	
	lim.glm<-glm(max.prey.size~length,data=list(max.prey.size=(max.prey.size),length=(as.numeric(names(max.prey.size)))))
	lim.glm2<-glm(max.prey.size~0+length,data=list(max.prey.size=(max.prey.size[max.prey.size>0]),
			length=(as.numeric(names(max.prey.size[max.prey.size>0])))),family=Gamma(link="log"))
	sp.dat$lim.glm.a[predd]<-(as.numeric(coef(lim.glm)[1]))
	sp.dat$lim.glm.b[predd]<-(as.numeric(coef(lim.glm)[2]))
	for(preyy in 1:2){
		
		pop<-tapply(subb$Wt,subb$Length,length)
		subb.l$ppratio<-subb.l$PREY_SIZE_CM/subb.l$PRED_LEN
		subb.l$logppratio<-log(subb.l$PREY_SIZE_CM/subb.l$PRED_LEN)
		###########################
		## Find max prey size
		###########################
		#max.prey.size<-tapply(subb.l$PREY_SIZE_CM,subb.l$PRED_LEN,max)		
		#max.prey.size<-tapply(subb.l$PREY_SIZE_CM[subb.l$ECOPATH_PREY==sp.dat$spp.l[preyy]],subb.l$PRED_LEN[subb.l$ECOPATH_PREY==sp.dat$spp.l[preyy]],max)
		#max.prey.size2<-tapply(subb.l$PREY_SIZE_CM[subb.l$ECOPATH_PREY==sp.dat$spp.l[preyy]],round(subb.l$ppratio[subb.l$ECOPATH_PREY==sp.dat$spp.l[preyy]],2),max)
		#max.prey.size<-max.prey.size[-1]
		plot(as.numeric(names(max.prey.size)),max.prey.size,col="red",pch=16)
		points(subb.l$PRED_LEN[subb.l$ECOPATH_PREY==sp.dat$spp.l[preyy]],subb.l$PREY_SIZE_CM[subb.l$ECOPATH_PREY==sp.dat$spp.l[preyy]])		
		points(as.numeric(names(max.prey.size)),max.prey.size,col="red",pch=16)
		#max.prey.size<-max.prey.size[max.prey.size>0]
		#lim.glm<-glm(max.prey.size~length,data=list(max.prey.size=(max.prey.size),length=(as.numeric(names(max.prey.size)))))
		lim.glm2<-glm(max.prey.size~0+length,data=list(max.prey.size=(max.prey.size[max.prey.size>0]),
				length=(as.numeric(names(max.prey.size[max.prey.size>0])))),family=Gamma(link="log"))
		#sp.dat$lim.glm.a[predd,preyy]<-exp(as.numeric(coef(lim.glm)[1]))
		#sp.dat$lim.glm.b[predd,preyy]<-(as.numeric(coef(lim.glm)[2]))
		lim<-predict.glm(lim.glm,type="response")
		gamma.shape(lim.glm2)
		lim2<-predict.glm(lim.glm2,type="response")
		lines(as.numeric(names(max.prey.size)),lim,col="red",lwd=2)		
		lines(as.numeric(names(max.prey.size[max.prey.size>0])),lim2,col="blue",lwd=2)		
		LL=1:180
		switch1<-rep(1,length(LL))
		switch2<-switch1
		prey.prop<-pop/sum(pop)
		LL.pred<-1:180
		LL.prey<-as.numeric(names(prey.prop))
		Kprop<-t(exp(sp.dat$a.glm[predd]+(LL.pred*sp.dat$b.glm[predd])))
		Kprop[Kprop>maxK[predd,preyy]]<-maxK[predd,preyy]
		if(preyy==1){Kprop[Kprop>maxK[predd,preyy]]<-0}
		Upaij.hat<-matrix(0,length(LL.prey),length(LL.pred))
		rownames(Upaij.hat)<-LL.prey
		colnames(Upaij.hat)<-LL.pred
		for(i in 1:length(LL.pred)){
			tlim<-(sp.dat$lim.glm.a[predd]+(LL.pred[i]*sp.dat$lim.glm.b[predd]))
			prey.prop.use<-prey.prop
			switch1<-rep(1,length(LL.prey))
			switch1[LL.prey>=tlim]<-exp(-(LL.prey[LL.prey>=tlim]-tlim)/10)
			prey.use<-prey.prop*switch1/sum(prey.prop*switch1)
			Upaij.hat[,i]<-prey.use*Kprop[i]
		}
		eval(parse(text=paste("Upaij.hat.",predd,preyy,"<-Upaij.hat",sep="")))	
	}
	rm(subb)
		rm(subb.l)
}
stomLength<-1:140
npred<-3
nsize<-length(stomLength)
StomWt<-matrix(0,npred,nsize)
StomWt[1,match(S.1$Length,stomLength)]<-S.1$S
StomWt[2,match(S.2$Length,stomLength)]<-S.2$S
StomWt[3,match(S.3$Length,stomLength)]<-S.3$S


#########################################################
## 3. Calculate ration for each pred - for now using livingston method where annual ration = 24*mean(S)*R and R = a*e^0.115*T, where S is stomach weight as %BW
#########################################################
Temp.data$R<-0.0143*exp(0.115*Temp.data$EBS_bottomT)*24
Temp.data$fT1<-bioenergetics(par=model.parms[1,],W=S.1$Wt[1],TempC<-Temp.data$EBS_bottomT,P=1,Eprey=3000,Efish=3000)$fTc
Temp.data$fT2<-bioenergetics(par=model.parms[4,],W=S.2$Wt[1],TempC<-Temp.data$EBS_bottomT,P=1,Eprey=3000,Efish=3000)$fTc
Temp.data$fT3<-bioenergetics(par=model.parms[5,],W=S.2$Wt[1],TempC<-Temp.data$EBS_bottomT,P=1,Eprey=3000,Efish=3000)$fTc

test<-bioenergetics(par=model.parms[1,],W=S.1$Wt,TempC<-Temp.data$EBS_bottomT[y],P=1,Eprey=3000,Efish=3000)
# ADMB code --> R[yr]*S[pred.l,prey.length]
C.liv.1<-Temp.data$R%*%t(S.1$S)
C.liv.2<-Temp.data$R%*%t(S.2$S)
C.liv.3<-Temp.data$R%*%t(S.3$S)

C.bioen.1<-Temp.data$fT1%*%t(S.1$C.hat.noFT)
C.bioen.2<-Temp.data$fT2%*%t(S.2$C.hat.noFT)
C.bioen.3<-Temp.data$fT3%*%t(S.3$C.hat.noFT)

#########################################################
## 4. Predict size at age using multinomial logistic regression regression (pcod and arrowtooth only) - need to specify for arrowooth...
#########################################################
graphics.off()
par(mfrow=c(1,3))
########### Pollock
dat<-plk_WA
dat2<-plk_WAT
ages<-sort(unique(dat$age))
#m.p<-glm(age~length,data=list(age=dat$age,length=dat$L),family=poisson(link = "log"))
m.p<-polr(factor(age)~length,data=list(age=dat$age,length=dat$L),method="logistic")  # ordered probit model - proportional odds logistic regression
m.pT<-polr(factor(age)~length*T,data=list(age=dat2$age,length=dat2$L,T=dat2$T),method="logistic")  # ordered probit model - proportional odds logistic regression
t.m.p<-data.frame(summary(m.p)[1])
plot(dat$age,dat$L, pch=16,cex=.8, main="pollock")
LatA.plk<-cbind(t.m.p[,1]/t.m.p[1,1],ages)
LatA.plk[,1]<-LatA.plk[,1]/10  # convert to cm

pollock.l$pred.age<-predict(m.p,newdata=data.frame(length=pollock.l$PRED_LEN*10))
pollock$pred.age<-as.numeric(predict(m.p,newdata=data.frame(length=pollock$Length*10)))
ages.pollock<-ages
tt<-tapply(pollock$pred.age,pollock$Length,mean.na)
points(pollock$pred.ageT,pollock$Length*10,col=colors()[111],pch=16)
points(pollock$pred.age,pollock$Length*10,col="red",pch=16)
points(t.m.p[-1,1]/t.m.p[1,1],ages[-1],type="b",pch=16,col="red",lwd=1,main="pollock")
lines(tt,as.numeric(names(tt))*10,col="red",lwd=2)
pollock$pred.ageT<-as.numeric(predict(m.pT,newdata=data.frame(length=pollock$Length*10,T=pollock$GearTemp)))
mp.plk<-m.p
LL<-seq(0.0001,1000,1)
LA.plk<-data.frame(Length=LL,age=as.numeric(predict(mp.plk,newdata=data.frame(length=LL))))
### MISC
#m.pw<-polr(factor(age)~weight,data=list(age=dat$age,weight=dat$W),method="logistic")  # ordered probit model - proportional odds logistic regression
#m.pwT<-polr(factor(age)~weight*T,data=list(age=dat2$age,weight=dat2$W,T=dat2$T),method="logistic")  # ordered probit model - proportional odds logistic regression
#t.m.pw<-data.frame(summary(m.pw)[1])
#plot(dat$age,dat$W, pch=16,cex=.8, main="pollock")
#AatW.plk<-cbind(ages,t.m.pw[,1]/t.m.pw[1,1])
#ages.pollock<-ages
#points(ages[-1],t.m.pw[-1,1]/t.m.pw[1,1],type="b",pch=16,col="red",lwd=1,main="pollock")
##VONB
#outfile<-"/Users/kkari/Documents/science/Projects/MSM/KerimVonB/atf/s_vonb.dat"
#results<-reptoRlist("pollock/s_vonb.rep")
#results.plk<-results
##plot(results$age,exp(results$size)/1000,ylab="Weight (Kg)",main="Pollock",xlab="Age (years)",pch=16,cex=.8)
#points(results$age,exp(results$pred_size),ylab="Weight (Kg)",xlab="Age (years)",pch=16,cex=1.5,col=colors()[111])
#tt<-tapply(dat$W,dat$age,mean.na)
#lines(as.numeric(names(tt)),tt,col="red",lwd=2)


########### PCOD
dat<-pcod_WA
dat2<-pcod_WAT
ages<-sort(unique(dat$age))
#m.p<-glm(age~length,data=list(age=dat$age,length=dat$L),family=poisson(link = "log"))
m.p<-polr(factor(age)~length,data=list(age=dat$age,length=dat$L),method="logistic")  # ordered probit model - proportional odds logistic regression
m.pT<-polr(factor(age)~length*T,data=list(age=dat2$age,length=dat2$L,T=dat2$T),method="logistic")  # ordered probit model - proportional odds logistic regression
t.m.p<-data.frame(summary(m.p)[1])
plot(dat$age,dat$L, pch=16,cex=.8, main="pcod")
LatA.pcod<-cbind(t.m.p[,1]/t.m.p[1,1],ages)
LatA.pcod[,1]<-LatA.pcod[,1]/10  # convert to cm
pcod.l$pred.age<-predict(m.p,newdata=data.frame(length=pcod.l$PRED_LEN*10))
pcod$pred.age<-as.numeric(predict(m.p,newdata=data.frame(length=pcod$Length*10)))
ages.pcod<-ages
tt<-tapply(pcod$pred.age,pcod$Length,mean.na)
points(pcod$pred.ageT,pcod$Length*10,col=colors()[111],pch=16)
points(pcod$pred.age,pcod$Length*10,col="red",pch=16)
points(t.m.p[-1,1]/t.m.p[1,1],ages[-1],type="b",pch=16,col="red",lwd=1,main="pcod")
lines(tt,as.numeric(names(tt))*10,col="red",lwd=2)
pcod$pred.ageT<-as.numeric(predict(m.pT,newdata=data.frame(length=pcod$Length*10,T=pcod$GearTemp)))

mp.pcod<-m.p
LL<-seq(0.0001,1000,1)
LA.pcod<-data.frame(Length=LL,age=as.numeric(predict(mp.pcod,newdata=data.frame(length=LL))))
########### ATF
dat<-atf_WA
dat2<-atf_WAT
ages<-sort(unique(dat$age))
#m.p<-glm(age~length,data=list(age=dat$age,length=dat$L),family=poisson(link = "log"))
m.p<-polr(factor(age)~length,data=list(age=dat$age,length=dat$L),method="logistic")  # ordered probit model - proportional odds logistic regression
m.pT<-polr(factor(age)~length*T,data=list(age=dat2$age,length=dat2$L,T=dat2$T),method="logistic")  # ordered probit model - proportional odds logistic regression
t.m.p<-data.frame(summary(m.p)[1])
plot(dat$age,dat$L, pch=16,cex=.8, main="arrowtooth")
LatA.arrowtooth<-cbind(t.m.p[,1]/t.m.p[1,1],ages)
LatA.arrowtooth[,1]<-LatA.arrowtooth[,1]/10  # convert to cm
arrowtooth.l$pred.age<-predict(m.p,newdata=data.frame(length=arrowtooth.l$PRED_LEN*10))
arrowtooth$pred.age<-as.numeric(predict(m.p,newdata=data.frame(length=arrowtooth$Length*10)))
ages.arrowtooth<-ages
tt<-tapply(arrowtooth$pred.age,arrowtooth$Length,mean.na)
points(arrowtooth$pred.ageT,arrowtooth$Length*10,col=colors()[111],pch=16)
points(arrowtooth$pred.age,arrowtooth$Length*10,col="red",pch=16)
points(t.m.p[-1,1]/t.m.p[1,1],ages[-1],type="b",pch=16,col="red",lwd=1,main="arrowtooth")
lines(tt,as.numeric(names(tt))*10,col="red",lwd=2)
arrowtooth$pred.ageT<-as.numeric(predict(m.pT,newdata=data.frame(length=arrowtooth$Length*10,T=arrowtooth$GearTemp)))
mp.arrowtooth<-m.p
LL<-seq(0.0001,1000,1)
LA.arrowtooth<-data.frame(Length=LL,age=as.numeric(predict(mp.arrowtooth,newdata=data.frame(length=LL))))
plot(LA.plk[,2],LA.plk[,1],type="l",lwd=2)
points(LA.pcod[,2],LA.pcod[,1],type="l",col="red",lwd=2)
points(LA.arrowtooth[,2],LA.arrowtooth[,1],type="l",col="blue",lwd=2)


################################################################
## Figure of K rations and fits
################################################################
dat1<-rep(0,6)
for(sp in 1:3){
	sub.dat<-prey[prey$Species==sp.dat$spp.p.pred[sp],]
	tmp<-data.frame(cbind(GearTemp=sub.dat$GearTemp,Length=sub.dat$Length,Total=sub.dat$TotWt,pollock=sub.dat$WALLEYE.POLLOCK_Wt,pcod=sub.dat$PACIFIC.COD_Wt))
	tmp$pred=rep(sp.dat$spp.p.pred[sp],length(tmp[,1]))
	dat1<-rbind(dat1,tmp)
}
dat1<-dat1[-1,]
dat1$Total[dat1$Total==0]<-NA
dat1$plk.prop=dat1$pollock/dat1$Total
dat1$pcod.prop=dat1$pcod/dat1$Total
dat1$other.prop=1-apply(cbind(dat1$plk.prop,dat1$pcod.prop),1,sum.na)
dat1.all<-dat1

dev.new(width=6.5,height=5)
	layout(rbind(c(1,0,3,0,5),
				c(2,0,4,0,6))
		,widths=c(1,0.01,1,0.01,1),
		heights=rep(1,2))
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(5,5,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(6)
## derive parameters for the gamma function for proportion of each prey item in each pred diet
maxK<-matrix(1,n.spp,n.spp)
maxKL<-maxK
predd<-1;preyy<-1

lbins<-list(
	plk<-c(14,24,32,39,43,48,51,54,57,59,60,61),
	pcod<-c(9,12,15,18,21,24,27,30,33,36,39,42,45,50,55,60,65,70,75,80,85,90,95,100,105),
	atf<-c(10,16,18,20,22,24,26,28,30,32,34,36,38,40,43,46,49,52,55,58,61,64,67,70,75)) # 95
nlengths<-c(12,25,25)
Kuse<-list(
	plk<-matrix(0,2,nlengths[1]),
	pcod<-matrix(0,2,nlengths[2]),
	atf<-matrix(0,2,nlengths[3]))
mm<-1:130	
Kuse<-list(
	plk<-matrix(0,2,length(mm)),
	pcod<-matrix(0,2,length(mm)),
	atf<-matrix(0,2,length(mm)))
	
predd<-1
preyy<-1

for (predd in 1:n.spp){
	for (preyy in 1:(n.spp-1)){
		cc<-preyy+6
		dat1<-dat1.all[dat1.all$pred==sp.dat$spp.p.pred[predd],]
		dat1.2<-dat1.all[dat1.all$pred==sp.dat$spp.p.pred[predd],]
		dat2.2<-na.omit(data.frame(cbind(Length=dat1$Length,tot=dat1$Total,plk=dat1[,cc])))
		#dat1<-na.omit(dat1[dat1$plk.prop>0,])
		dat2<-na.omit(data.frame(cbind(Length=dat1$Length,prop=dat1[,cc])))
		dat2$prop<-dat2$prop+1e-2
		#dat2<-dat2[dat2$prop>0,]
		mdat<-tapply(dat2$prop,dat2$Length,mean.na)-1e-2
		ldat<-tapply(dat2$prop,dat2$Length,length.na)
		ldat<-ldat/sum(ldat)
		llim<-c(as.numeric(names(ldat))[ldat<0.0001],max(dat2$Length))
		limm<-min(llim[llim>50])
		dat2old<-dat2
		dat2<-dat2[dat2$Length<=limm,]
		#points(as.numeric(names(mdat)),mdat,lwd=2,pch=16)
		#plot(density(dat2$prop+1e-2),xlim=c(0,max(mdat)*1.1))
		#gam.1<-gam(cbind(dat2.2$total-dat2.2$plk,dat2.2$plk)~s(dat2.2$Length,1),family=binomial(link="log"))
		
		model.a<-glm(prop~Length+I(Length^2)+I(Length^3)+I(Length^4)+I(Length^5)+I(Length^6),data=dat2,family=Gamma(link="log"))

		gam.1<-gam(prop~s(Length,1),family=Gamma(link="log"),data=dat2)
		gam.2<-gam(prop~s(dat2$Length,2),family=Gamma(link="log"),data=dat2)
		gam.3<-gam(prop~s(Length,3),family=Gamma(link="log"),data=dat2)
		gam.4<-gam(prop~s(Length,4),family=Gamma(link="log"),data=dat2)
		gam.5<-gam(prop~s(Length,5),family=Gamma(link="log"),data=dat2)
		gam.6<-gam(prop~s(Length,6),family=Gamma(link="log"),data=dat2)
		gam.7<-gam(prop~s(Length,7),family=Gamma(link="log"),data=dat2)
		gam.8<-gam(prop~s(Length,8),family=Gamma(link="log"),data=dat2)
				
		aa<-AIC(gam.1,gam.2,gam.3,gam.4,gam.5,gam.6,gam.7,gam.8)
		aa$num<-1:length(aa)
		aa$deltaAIC<-aa$AIC-min(aa$AIC)
		aa$rank<-rank(aa$deltaAIC)
		aa<-aa[order(aa$rank),]
		eval(parse(text=paste("gam.use<-",rownames(aa)[1],sep="")))
		eval(parse(text=paste("gam.use",predd,preyy,"<-gam.use",sep="")))
		eval(parse(text=paste("glm.use",predd,preyy,"<-model.a",sep="")))
		llengths<-seq(0,120,.5)
		#Kuse[[predd]][preyy,]<-exp(predict.gam(gam.use,newdata=list(Length=lbins[[predd]])))-1e-2
		Kuse[[predd]][preyy,]<-exp(predict.gam(gam.use,newdata=list(Length=mm)))-1e-2
		
		tt2<-(predict.gam(gam.use,se.fit = TRUE))
		tt2.m<-tapply(exp(tt2$fit),dat2$Length,mean.na)-1e-2
		tt2.u<-tapply(exp(tt2$fit+1.95*tt2$se.fit),dat2$Length,mean.na)-1e-2
		tt2.l<-tapply(exp(tt2$fit-1.95*tt2$se.fit),dat2$Length,mean.na)-1e-2
		tt4<-exp(predict.gam(gam.use,newdata=list(Length=1:130)))-1e-2
		xx<-names(tt2.m)
		nx<-length(xx)
		cexx<-.7
		plot(as.numeric(names(mdat)),mdat,lwd=2,pch=16,ylim=c(0,max(mdat)*1.1),col="white",axes=FALSE);axis(1);axis(2)
		eval(parse(text=paste("gam.use<-",rownames(aa)[1],sep="")))
		mtext(side=3,rownames(aa)[1],line=-1)
		
		points(dat2$Length,dat2$prop-1e-2,pch=16,col=colors()[240],cex=cexx)
		polygon(c(xx,xx[nx:1]),c(tt2.u,tt2.l[nx:1]),border=FALSE,col=colors()[200])
		points(as.numeric(names(mdat)),mdat,lwd=1,pch=16,cex=cexx,type="p")
		lines(names(tt2.m),tt2.m,col=colors()[24],lwd=2)
		lines(names(tt4),tt4,col=colors()[24],lwd=2)
	}		
}

Kuse[[1]][Kuse[[1]]<0]<-0
Kuse[[3]][Kuse[[3]]<0]<-0
Kuse[[2]][Kuse[[2]]<0]<-0



#########################################################
## 5. Remove zeros and output L at Age to .dat file for K parms
#########################################################
outfile<-"~/Documents/science/SVN/MSM/svn (trunk)/KirOct2011/K.dat"
# K regression coef for each sppXspp combo where K <-a*pred_size^b
cat("# K_Coef for K size relationship ",file=outfile,append=FALSE,sep="\n")

cat(coef(glm.use11),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat(coef(glm.use12),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")

cat(coef(glm.use21),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat(coef(glm.use22),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")

cat(coef(glm.use31),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat(coef(glm.use32),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")


#########################################################
## 5. Remove zeros and output L at Age to .dat file
#########################################################
outfile<-"/Users/kkari/Documents/science/Projects/MSM/msm_data/stomdat/stomach1.dat"
outfile<-"~/Documents/science/SVN/MSM/svn (trunk)/KirOct2011/stomach1.dat"
updateFig<-1
outfile2<-"~/Documents/science/SVN/MSM/svn (trunk)/KirOct2011/MSM_figures/"

# K regression coef for each sppXspp combo where K <-a*pred_size^b
cat("# # Number of data points ",file=outfile,append=FALSE,sep="\n")
cat(length(datause[,1]),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# ages ",file=outfile,append=TRUE,sep="\n")
cat(as.numeric(LA.plk$age),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
cat(as.numeric(LA.pcod$age),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
cat(as.numeric(LA.arrowtooth$age),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Length ",file=outfile,append=TRUE,sep="\n")
cat(as.numeric(LA.plk$Length),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
cat(as.numeric(LA.pcod$Length),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
cat(as.numeric(LA.arrowtooth$Length),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")

cat("// Klengths: lengths for the Khat ",file=outfile,append=TRUE,sep="\n")
cat(mm,file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("// Khat: Proportion of total stomach by length ",file=outfile,append=TRUE,sep="\n")

for(predd in 1:3){
cat(Kuse[[predd]][1,],file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat(Kuse[[predd]][2,],file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
}
#########################################################
## 5. Remove zeros and output to .dat file
#########################################################
save(list=ls(),file="pseudoCode.RData")
setwd(savefile)

cat("// Tyrs: Temperature years ",file=outfile,append=app1,sep="\n")
cat(t(Temp.data$Year),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("// TempC: Temperature ",file=outfile,append=TRUE,sep="\n")
cat(t(Temp.data$EBS_bottomT),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("// R: R from Livingston ",file=outfile,append=TRUE,sep="\n")
cat(t(Temp.data$R),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("// fT: Temperature function from bioenergetics models ",file=outfile,append=TRUE,sep="\n")
cat(t(Temp.data$fT1),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat(t(Temp.data$fT2),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat(t(Temp.data$fT3),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("// nTyrs: number of Temperature years ",file=outfile,append=TRUE,sep="\n")
cat(length(Temp.data$Year),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("// nTyrs: number of Temperature years ",file=outfile,append=TRUE,sep="\n")
cat(length(Temp.data$Year),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")


##################################################################################################################
##################################################################################################################
##################################################################################################################


#########################################################
## Fig 1. K as a function of predator size
#########################################################
graphics.off()
dev.new(width=8.5,height=6)
	layout(rbind(c(1,0,3,0,5),
				c(2,0,4,0,6))
		,widths=c(1,0.01,1,0.01,1),
		heights=rep(1,2))
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(5,5,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(6)
predd<-1
preyy<-1
usegam<-0
for (predd in 1:n.spp){
	for (preyy in 1:(n.spp-1)){
		cc<-preyy+6
		dat1<-dat1.all[dat1.all$pred==sp.dat$spp.p.pred[predd],]
		dat1.2<-dat1.all[dat1.all$pred==sp.dat$spp.p.pred[predd],]
		dat2.2<-na.omit(data.frame(cbind(Length=dat1$Length,tot=dat1$Total,plk=dat1[,cc])))
		dat2<-na.omit(data.frame(cbind(Length=dat1$Length,prop=dat1[,cc])))
		dat2$prop<-dat2$prop+1e-2
		mdat<-tapply(dat2$prop,dat2$Length,mean.na)-1e-2
		ldat<-tapply(dat2$prop,dat2$Length,length.na)
		ldat<-ldat/sum(ldat)
		llim<-c(as.numeric(names(ldat))[ldat<0.0001],max(dat2$Length))
		limm<-min(llim[llim>50])
		dat2old<-dat2
		dat2<-dat2[dat2$Length<=limm,]
		if(usegam==1){
			eval(parse(text=paste("gam.use<-gam.use",predd,preyy,sep="")))
		}else{
			eval(parse(text=paste("gam.use<-glm.use",predd,preyy,sep="")))
		}
		tt2<-(predict.gam(gam.use,se.fit = TRUE))
		tt2.m<-tapply(exp(tt2$fit),dat2$Length,mean.na)-1e-2
		tt2.u<-tapply(exp(tt2$fit+1.95*tt2$se.fit),dat2$Length,mean.na)-1e-2
		tt2.l<-tapply(exp(tt2$fit-1.95*tt2$se.fit),dat2$Length,mean.na)-1e-2
		tt4<-exp(predict.gam(gam.use,newdata=list(Length=1:130)))-1e-2
		xx<-names(tt2.m)
		nx<-length(xx)
		cexx<-.7
		plot(as.numeric(names(mdat)),mdat,lwd=2,pch=16,ylim=c(0,max(mdat)*1.1),col="white",axes=FALSE)
		axis(1,at=seq(0,120,20));axis(2)
		points(dat2$Length,dat2$prop-1e-2,pch=16,col=colors()[240],cex=cexx)
		polygon(c(xx,xx[nx:1]),c(tt2.u,tt2.l[nx:1]),border=FALSE,col=colors()[200])
		points(as.numeric(names(mdat)),mdat,lwd=1,pch=16,cex=cexx,type="p")
		lines(names(tt2.m),tt2.m,col=colors()[24],lwd=2)
		lines(names(tt4),tt4,col=colors()[24],lwd=2)
		mtext(side=3,paste(sp.dat$spp[predd],"-->",sp.dat$spp[preyy]),line=-2)		

	}		
}
		mtext(side=1,"Predator size (cm)",font=2,line=2,outer=TRUE)		
		mtext(side=2,"Stomach proportion by weight (K)",font=2,line=2,outer=TRUE)		
if(updateFig==1){dev.copy2pdf(device = quartz, file =paste(outfile2,"Kplot.pdf",sep=""))}

#########################################################
## Fig 2. g/g as a function of predator size
#########################################################

vonb_fun<-function(par,data){
	rm(size)
	rm(age)
	rm(data)
	data<-data
	age<-data$age
 	size<-data$L
	Winf<-1.1*max((size)) # // set linf to 1.1 times the longest observed length
  	
  	K<-par[1]
  	t0<-par[2]
  	dexp<-par[3]
  	sigma<-exp(par[4])
  	
  	pred_size<-((Winf)+(1/(1-dexp))*log(1-exp(-K*(1-dexp)*(age-t0))))
    H<-K*(Winf^(1-dexp))
    likeH<-H
    nobs<-length(age)
    f<-sum((log(size)-log(pred_size))^2)
    # f<-norm2(pred_size-size);
    f<-nobs*log(sigma) + f/(2.0*sigma*sigma)
    #f = regression(size,pred_size);
  	return(f)
}

graphics.off()
Pval<-c(2,1,3)
X11(width=8.5,height=6)
	layout(rbind(c(1,0),
				c(2,0),
				c(3,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(3)
col1<-colors()[351]
col2<-colors()[136]
col3<-colors()[496]
col4<-colors()[638]

predd<-2
preyy<-1
spnames<-c("plk","pcod","atf")
for(predd in 1:3){
	rm(data)
	rm(vonb)
	eval(parse(text=paste("data<-",spnames[predd],"_WA",sep="")))
	vonb<-optim(vonb_fun,par=c(0.5,-0.25,0.8,log(1.0)),data=data)
	eval(parse(text=paste("vonb",predd,"<-vonb",sep="")))
}

for (predd in 1:n.spp){
		cc<-preyy+5
		dat1<-dat1.all[dat1.all$pred==sp.dat$spp.p.pred[predd],]
		dat1$Wt<-sp.dat$LW.a[predd]*dat1$L^sp.dat$LW.b[predd]
		eval(parse(text=paste("LWdat<-LW_",spnames[predd],sep="")))
		#plot(LWdat$L,LWdat$W)
		#points(dat1$L,dat1$Wt,col="red")
		dat1$specC<-dat1$Total/dat1$Wt
		dat1$digC<-dat1$Total*24*0.0134*exp(.115*dat1$GearTemp)/dat1$Wt
		dat1$oldMSM<-dat1$Total*24*0.0134*exp(.115*3)/dat1$Wt
		parms<-model.parms.sub[predd,]
		C<-Consum_fun(parms=parms,W=dat1$Wt,TempC=dat1$GearTemp)*Pval[predd]
		dat1$BioenC<-C
		#Kerim's code
		
		# Indigestable stuff is 40%
		 AAcons <- 0.6
		eval(parse(text=paste("vonb<-vonb",predd,sep="")))
		size<-data$L
		Winf<-1.1*max((size)) # // set linf to 1.1 times the longest observed length	
		K<-vonb$par[1]
		t0<-vonb$par[2]
		dexp<-vonb$par[3]
		sigma<-vonb$par[4]
		H<-K*(Winf^(1-dexp))		
		AAcons <- 0.6
		lenlist<- dat1$Length
		wlist<- dat1$Wt
		dailycons<-(H*(wlist^dexp)/AAcons)/365
		dat1$VonB<-dailycons/wlist
  		dat2<-na.omit(dat1)		
		
		mdat<-tapply(dat2$specC,dat2$Length,mean.na)
		mdat2<-tapply(dat2$digC,dat2$Length,mean.na)
		mdat.oldmsm<-tapply(dat2$oldMSM,dat2$Length,mean.na)
		mdat.bion<-tapply(dat2$BioenC,dat2$Length,mean.na)
		mdat.vonb<-tapply(dat2$VonB,dat2$Length,mean.na)

		se.dat<-tapply(dat2$specC,dat2$Length,se.na)
		mdat.u<-mdat+1.95*se.dat
		mdat.l<-mdat-1.95*se.dat
		ldat<-tapply(dat2$specC,dat2$Length,length.na)
		ldat<-ldat/sum(ldat)
		
		llim<-c(as.numeric(names(ldat))[ldat<0.0001],max(dat2$Length))
		limm<-min(llim[llim>50])
		dat2old<-dat2
		dat2<-dat2[dat2$Length<=limm,]
		cexx<-.9
		plot(names(mdat),mdat,pch=16,type="b",cex=cexx,ylab="",xlab="",main=sp.dat$spp[predd],axes=FALSE,ylim=c(0,1.1*max(mdat)));axis(1);axis(2)
		points(dat2$Length,dat2$specC,pch=16,col=col1,cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],ylim=c(0,.06),line=-1,,axes=FALSE);axis(1);axis(2)
		xx<-as.numeric(names(mdat))
	
		ns<-length(mdat)
		xx<-c(xx,xx[ns:1])
		yy<-c(mdat.u,mdat.l[ns:1])
		polydat<-na.omit(data.frame(xx,yy))
		polygon(polydat[,1],polydat[,2],border=FALSE,col=colors()[220])
		points(names(mdat),mdat,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],axes=FALSE,lwd=2);axis(1);axis(2)
		points(names(mdat2),mdat2,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col3)
		points(names(mdat.oldmsm),mdat.oldmsm,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col4)
		points(names(mdat.bion),mdat.bion,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col2)
		points(names(mdat.vonb),mdat.vonb,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col2,lty=2)

		lines(1:40,rep(.41/100,length(1:40)),col=col3,lwd=2,lty=1)
		lines(40:130,rep(.46/100,length(40:130)),col=col3,lwd=2,lty=1)
		#lines(S.1$Length,S.1$C,pch=16,col=colors()[300],cex=1,lty=1,lwd=2)
		if(predd==1){
		legend(45,.07,c("mean data (Stom.*digest)","VonB","bioenergetics","livingston","old MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
		}	
	}

		mtext(side=1,"Predator size (cm)",font=2,line=2,outer=TRUE)		
		mtext(side=2,"Specific daily consumption rate (g/g)",font=2,line=2,outer=TRUE)		
if(updateFig==1){dev.copy2pdf(device = quartz, file =paste(outfile2,"Cbyday.pdf",sep=""))}



X11(width=8.5,height=6)
	layout(rbind(c(1,0),
				c(2,0),
				c(3,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(3)
col1<-colors()[351]
col2<-colors()[136]
col3<-colors()[496]
col4<-colors()[638]

predd<-1
preyy<-1
spnames<-c("plk","pcod","atf")
model.parms.sub<-model.parms[c(1,4,5),]
vonbK<-list(Winf=21448.1, K=0.57325, t0=-0.495495, dexp=0.742355, H=7.48695)

for (predd in 1:n.spp){
		cc<-preyy+5
		dat1<-dat1.all[dat1.all$pred==sp.dat$spp.p.pred[predd],]
		dat1$Wt<-sp.dat$LW.a[predd]*dat1$L^sp.dat$LW.b[predd]
		eval(parse(text=paste("LWdat<-LW_",spnames[predd],sep="")))
		#plot(LWdat$L,LWdat$W)
		#points(dat1$L,dat1$Wt,col="red")
		dat1$specC<-dat1$Total
		dat1$digC<-dat1$Total*24*0.0134*exp(.115*dat1$GearTemp)
		dat1$oldMSM<-dat1$Total*24*0.0134*exp(.115*3)
		parms<-model.parms.sub[predd,]
		C<-Consum_fun(parms=parms,W=dat1$Wt,TempC=dat1$GearTemp)*Pval[predd]
		dat1$BioenC<-C*dat1$Wt
		
		
		#Kerim's code
		
		# Indigestable stuff is 40%
		 AAcons <- 0.6
		eval(parse(text=paste("vonb<-vonb",predd,sep="")))
		size<-data$L
		Winf<-1.1*max((size)) # // set linf to 1.1 times the longest observed length	
		K<-vonb$par[1]
		t0<-vonb$par[2]
		dexp<-vonb$par[3]
		sigma<-vonb$par[4]
		H<-K*(Winf^(1-dexp))		
		AAcons <- 0.6
		lenlist<- dat1$Length
		wlist<- dat1$Wt
		dailycons<-(H*(wlist^dexp)/AAcons)
		dat1$VonB<-dailycons
  		dat2<-na.omit(dat1)		
		
		
		if(predd==2){
			size<-data$L
			Winf<-1.1*max((size)) # // set linf to 1.1 times the longest observed length	
			K<-vonbK$K
			t0<-vonbK$t0
			dexp<-vonbK$dexp
			H<-vonbK$H
			Winf<-vonbK$Winf
			H<-K*(Winf^(1-dexp))		
			AAcons <- 0.6
			lenlist<- dat1$Length
			wlist<- dat1$Wt
			dailycons<-(H*(wlist^dexp)/AAcons)
			dat1$vonbKerim<-dailycons
		}
		dat2<-na.omit(dat1)	
		mdat<-tapply(dat2$specC,dat2$Length,mean.na)
		mdat2<-tapply(dat2$digC,dat2$Length,mean.na)
		mdat.oldmsm<-tapply(dat2$oldMSM,dat2$Length,mean.na)
		mdat.bion<-tapply(dat2$BioenC,dat2$Length,mean.na)
		mdat.vonb<-tapply(dat2$VonB,dat2$Length,mean.na)
		if(predd==2){
			mdat.vonbK<-tapply(dat2$vonbKerim,dat2$Length,mean.na)
		}

		se.dat<-tapply(dat2$specC,dat2$Length,se.na)
		mdat.u<-mdat+1.95*se.dat
		mdat.l<-mdat-1.95*se.dat
		ldat<-tapply(dat2$specC,dat2$Length,length.na)
		ldat<-ldat/sum(ldat)
		
		llim<-c(as.numeric(names(ldat))[ldat<0.0001],max(dat2$Length))
		limm<-min(llim[llim>50])
		dat2old<-dat2
		dat2<-dat2[dat2$Length<=limm,]
		cexx<-.9
		plot(names(mdat),mdat,pch=16,type="b",cex=cexx,ylab="",xlab="",main=sp.dat$spp[predd],axes=FALSE,ylim=c(0,1.1*max(mdat)));axis(1);axis(2)
		points(dat2$Length,dat2$specC,pch=16,col=col1,cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],ylim=c(0,.06),line=-1,,axes=FALSE);axis(1);axis(2)
		xx<-as.numeric(names(mdat))
	
		ns<-length(mdat)
		xx<-c(xx,xx[ns:1])
		yy<-c(mdat.u,mdat.l[ns:1])
		polydat<-na.omit(data.frame(xx,yy))
		polygon(polydat[,1],polydat[,2],border=FALSE,col=colors()[220])
		points(names(mdat),mdat,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],axes=FALSE,lwd=2);axis(1);axis(2)
		points(names(mdat2),mdat2,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col3)
		points(names(mdat.oldmsm),mdat.oldmsm,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col4)
		points(names(mdat.bion),mdat.bion,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col2)
		points(names(mdat.vonb),mdat.vonb/365,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col2,lty=2)
		points(names(mdat.vonbK),mdat.vonbK/365,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col="blue",lty=2)

		
		if(predd==1){
		lines(1:40,sp.dat$LW.a[predd]*(1:40)^sp.dat$LW.b[predd]*rep(.41/100,length(1:40)),col=col3,lwd=2,lty=1)
		lines(40:130,sp.dat$LW.a[predd]*(40:130)^sp.dat$LW.b[predd]*rep(.46/100,length(40:130)),col=col3,lwd=2,lty=1)
		}
		if(predd==2){
				lines(1:40,sp.dat$LW.a[predd]*(1:40)^sp.dat$LW.b[predd]*rep(.62/100,length(1:40)),col=col3,lwd=2,lty=1)
			lines(40:130,sp.dat$LW.a[predd]*(40:130)^sp.dat$LW.b[predd]*rep(1.33/100,length(40:130)),col=col3,lwd=2,lty=1)
		}
		#lines(S.1$Length,S.1$C,pch=16,col=colors()[300],cex=1,lty=1,lwd=2)
		if(predd==1){
		legend(45,.07,c("mean data (Stom.*digest)","VonB","bioenergetics","livingston","old MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
		}	
	}

		mtext(side=1,"Predator size (cm)",font=2,line=2,outer=TRUE)		
		mtext(side=2,"Daily ration (g/d)",font=2,line=2,outer=TRUE)		
if(updateFig==1){dev.copy2pdf(device = quartz, file =paste(outfile2,"dailyRation.pdf",sep=""))}




X11(width=8.5,height=6)
	layout(rbind(c(1,0),
				c(2,0),
				c(3,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(3)
col1<-colors()[351]
col2<-colors()[136]
col3<-colors()[496]
col4<-colors()[638]

predd<-1
preyy<-1
spnames<-c("plk","pcod","atf")
model.parms.sub<-model.parms[c(1,4,5),]
vonbK<-list(Winf=21448.1, K=0.57325, t0=-0.495495, dexp=0.742355, H=7.48695)

for (predd in 1:n.spp){
		cc<-preyy+5
		dat1<-dat1.all[dat1.all$pred==sp.dat$spp.p.pred[predd],]
		dat1$Wt<-sp.dat$LW.a[predd]*dat1$L^sp.dat$LW.b[predd]
		eval(parse(text=paste("LWdat<-LW_",spnames[predd],sep="")))
		#plot(LWdat$L,LWdat$W)
		#points(dat1$L,dat1$Wt,col="red")
		dat1$specC<-dat1$Total*91.5
		dat1$digC<-dat1$Total*24*0.0134*exp(.115*dat1$GearTemp)*91.5
		dat1$oldMSM<-dat1$Total*24*0.0134*exp(.115*3)*91.5
		parms<-model.parms.sub[predd,]
		C<-Consum_fun(parms=parms,W=dat1$Wt,TempC=dat1$GearTemp)*Pval[predd]
		dat1$BioenC<-C*dat1$Wt*91.5
		
		
		#Kerim's code
		
		# Indigestable stuff is 40%
		 AAcons <- 0.6
		eval(parse(text=paste("vonb<-vonb",predd,sep="")))
		size<-data$L
		Winf<-1.1*max((size)) # // set linf to 1.1 times the longest observed length	
		K<-vonb$par[1]
		t0<-vonb$par[2]
		dexp<-vonb$par[3]
		sigma<-vonb$par[4]
		H<-K*(Winf^(1-dexp))		
		AAcons <- 0.6
		lenlist<- dat1$Length
		wlist<- dat1$Wt
		dailycons<-(H*(wlist^dexp)/AAcons)
		dat1$VonB<-dailycons
  		dat2<-na.omit(dat1)		
		
		
		if(predd==2){
			size<-data$L
			Winf<-1.1*max((size)) # // set linf to 1.1 times the longest observed length	
			K<-vonbK$K
			t0<-vonbK$t0
			dexp<-vonbK$dexp
			H<-vonbK$H
			Winf<-vonbK$Winf
			H<-K*(Winf^(1-dexp))		
			AAcons <- 0.6
			lenlist<- dat1$Length
			wlist<- dat1$Wt
			dailycons<-(H*(wlist^dexp)/AAcons)
			dat1$vonbKerim<-dailycons
		}
		dat2<-na.omit(dat1)	
		mdat<-tapply(dat2$specC,dat2$Length,mean.na)
		mdat2<-tapply(dat2$digC,dat2$Length,mean.na)
		mdat.oldmsm<-tapply(dat2$oldMSM,dat2$Length,mean.na)
		mdat.bion<-tapply(dat2$BioenC,dat2$Length,mean.na)
		mdat.vonb<-tapply(dat2$VonB,dat2$Length,mean.na)
		if(predd==2){
			mdat.vonbK<-tapply(dat2$vonbKerim,dat2$Length,mean.na)
		}

		se.dat<-tapply(dat2$specC,dat2$Length,se.na)
		mdat.u<-mdat+1.95*se.dat
		mdat.l<-mdat-1.95*se.dat
		ldat<-tapply(dat2$specC,dat2$Length,length.na)
		ldat<-ldat/sum(ldat)
		
		llim<-c(as.numeric(names(ldat))[ldat<0.0001],max(dat2$Length))
		limm<-min(llim[llim>50])
		dat2old<-dat2
		dat2<-dat2[dat2$Length<=limm,]
		cexx<-.9
		plot(names(mdat),mdat,pch=16,type="b",cex=cexx,ylab="",xlab="",main=sp.dat$spp[predd],axes=FALSE,ylim=c(0,1.1*max(mdat)));axis(1);axis(2)
		points(dat2$Length,dat2$specC,pch=16,col=col1,cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],ylim=c(0,.06),line=-1,,axes=FALSE);axis(1);axis(2)
		xx<-as.numeric(names(mdat))
	
		ns<-length(mdat)
		xx<-c(xx,xx[ns:1])
		yy<-c(mdat.u,mdat.l[ns:1])
		polydat<-na.omit(data.frame(xx,yy))
		polygon(polydat[,1],polydat[,2],border=FALSE,col=colors()[220])
		points(names(mdat),mdat,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],axes=FALSE,lwd=2);axis(1);axis(2)
		points(names(mdat2),mdat2,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col3)
		points(names(mdat.oldmsm),mdat.oldmsm,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col4)
		points(names(mdat.bion),mdat.bion,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col2)
		points(names(mdat.vonb),mdat.vonb,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col=col2,lty=2)
		points(names(mdat.vonbK),mdat.vonbK,pch=16,type="l",cex=cexx,ylab="",xlab="",main=sp.dat$spp[1],lwd=2,col="blue",lty=2)

		
		if(predd==1){
		lines(1:40,91.5*sp.dat$LW.a[predd]*(1:40)^sp.dat$LW.b[predd]*rep(.41/100,length(1:40)),col=col3,lwd=2,lty=1)
		lines(40:130,91.5*sp.dat$LW.a[predd]*(40:130)^sp.dat$LW.b[predd]*rep(.46/100,length(40:130)),col=col3,lwd=2,lty=1)
		}
		if(predd==2){
				lines(1:40,91.5*sp.dat$LW.a[predd]*(1:40)^sp.dat$LW.b[predd]*rep(.62/100,length(1:40)),col=col3,lwd=2,lty=1)
			lines(40:130,91.5*sp.dat$LW.a[predd]*(40:130)^sp.dat$LW.b[predd]*rep(1.33/100,length(40:130)),col=col3,lwd=2,lty=1)
		}
		#lines(S.1$Length,S.1$C,pch=16,col=colors()[300],cex=1,lty=1,lwd=2)
		if(predd==1){
		legend(10,6000,c("mean Stom. Wt","VonB","bioenergetics","livingston","old MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
		}	
	}

		mtext(side=1,"Predator size (cm)",font=2,line=2,outer=TRUE)		
		mtext(side=2,"Annual ration (g/yr)",font=2,line=2,outer=TRUE)		
if(updateFig==1){dev.copy2pdf(device = quartz, file =paste(outfile2,"annualRation.pdf",sep=""))}
#########################################################
## 5. plot results
#########################################################
graphics.off()
# FIG 1

dev.new(width=4,height=6.5)
	layout(rbind(c(1,0),
				c(2,0),
				c(3,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(3)
subb<-S.1[S.1$S>0.005,]
plot(subb$Length,subb$S,pch=16,ylim=c(0,.1),ylab="",xlab="",xlim=c(0,80),main=sp.dat$spp[1],axes=FALSE);axis(1);axis(2)
points(subb$Length,subb$S.hat,type="l",col="red",lwd=2,ylim=c(0,.1))
points(subb$Length,subb$S.hat2,type="l",col="red",lwd=2,ylim=c(0,.1),lty=2)
lines(1:40,rep(.41/(100*24*0.0134*exp(0.115*3)),length(1:40)),col="blue",lwd=2)
lines(40:max(subb$Length),rep(.46/(100*24*0.0134*exp(0.115*3)),length(40:max(subb$Length))),col="blue",lwd=2)
legend(60,.08,c("data","bioenergetics","livingston"),col=c("black","red","blue"),lwd=2,box.lwd=0)
subb<-S.2[S.2$S>0.005,]
plot(subb$Length,subb$S,pch=16,ylim=c(0,.1),ylab="mean stomach wt (g/g)",xlab="Length",main=sp.dat$spp[2])
points(subb$Length,subb$S.hat,type="l",col="red",lwd=2,ylim=c(0,.1))
lines(1:55,rep(.62/(100*24*0.0134*exp(0.115*3)),length(1:55)),col="blue",lwd=2)
lines(55:max(subb$Length),rep(1.33/(100*24*0.0134*exp(0.115*3)),length(55:max(subb$Length))),col="blue",lwd=2)
subb<-S.3[S.3$S>0.005,]
plot(subb$Length,subb$S,pch=16,ylim=c(0,.1),ylab="mean stomach wt (g/g)",xlab="Length",main=sp.dat$spp[3])
points(subb$Length,subb$S.hat,type="l",col="red",lwd=2,ylim=c(0,.1))
lines(1:40,rep(.61/(100*24*0.0134*exp(0.115*3)),length(1:40)),col="blue",lwd=2)
lines(40:max(subb$Length),rep(.63/(100*24*0.0134*exp(0.115*3)),length(40:max(subb$Length))),col="blue",lwd=2)

mtext(bquote(bold("Length (cm)")),side=1,line=1,outer=TRUE)
mtext(bquote(bold("mean stomach wt (g/g)")),side=2,line=1,outer=TRUE)


# FIG 1 : plot of mean stom weight and MSM regression 
#graphics.off()

dev.new(width=4,height=6.5)
	layout(rbind(c(1,0),
				c(2,0),
				c(3,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(3)
col1<-colors()[351]
col2<-colors()[136]
col3<-colors()[496]
col4<-colors()[638]
col5<-colors()[300]
plot(pollock$Length,pollock$TotWt/pollock$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[1],ylim=c(0,.1),xlim=c(0,80),axes=FALSE);axis(1);axis(2)
points(S.1$Length,S.1$S,type="b",pch=16,col="black",cex=.7,lty=1,lwd=1)
lines(1:40,rep(.41/(100*24*0.0134*exp(0.115*3)),length(1:40)),col=col3,lwd=2,lty=1)
lines(40:130,rep(.46/(100*24*0.0134*exp(0.115*3)),length(40:130)),col=col3,lwd=2,lty=1)
legend(45,.08,c("mean data","livingston","MSM"),col=c("black",col3,col4),lwd=2,box.lwd=0,box.lty=0)
lines(S.1$Length,S.1$S.MSM,col=col4,lwd=2)
points(S.1$Length,S.1$S,type="b",pch=16,col="black",cex=.7,lty=1,lwd=1)
plot(pcod$Length,pcod$TotWt/pcod$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[2],ylim=c(0,.1),xlim=c(0,120),axes=FALSE);axis(1);axis(2)
points(S.2$Length,S.2$S,type="b",pch=16,col="black",cex=.7,lty=1,lwd=1)
lines(1:55,rep(.62/(100*24*0.0134*exp(0.115*3)),length(1:55)),col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/(100*24*0.0134*exp(0.115*3)),length(55:130)),col=col3,lwd=2,lty=1)
lines(S.2$Length,S.2$S.MSM,col=col4,lwd=2)
points(S.2$Length,S.2$S,type="b",pch=16,col="black",cex=.7,lty=1,lwd=1)
plot(arrowtooth$Length,arrowtooth$TotWt/arrowtooth$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[3],ylim=c(0,.1),xlim=c(0,80),axes=FALSE);axis(1);axis(2)
points(S.3$Length,S.3$S,type="b",pch=16,col="black",cex=.7,lty=1,lwd=1)
lines(1:40,rep(.61/(100*24*0.0134*exp(0.115*3)),length(1:40)),col=col3,lwd=2,lty=1)
lines(40:130,rep(.63/(100*24*0.0134*exp(0.115*3)),length(40:130)),col=col3,lwd=2,lty=1)
#lines(S.3$Length,S.3$C,pch=16,col="black",cex=1,lty=1,lwd=2)
S.a<-S.3
S.a<-S.a[S.a$Length<78,]
lines(S.3$Length,S.3$S.MSM,col=col4,lwd=2)
points(S.3$Length,S.3$S,type="b",pch=16,col="black",cex=.7,lty=1,lwd=1)
mtext(bquote(bold("Length (cm)")),side=1,line=1,outer=TRUE)
mtext(bquote(bold("Mean Stom proportion (g of Stom. weight/g of body weight)")),side=2,line=1,outer=TRUE)
dev.print(device=postscript, paste("fig1StomProp.eps",sep=""), onefile=FALSE, horizontal=FALSE)

# FIG 2 : plot of Ration and models (where data = mean S * digestion)
#graphics.off()

dev.new(width=4,height=6.5)
	layout(rbind(c(1,0),
				c(2,0),
				c(3,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(3)
col1<-colors()[351]
col2<-colors()[136]
col3<-colors()[496]
col4<-colors()[638]
plot(pollock$Length,pollock$TotWt/pollock$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[1],ylim=c(0,.06),line=-1,,axes=FALSE);axis(1);axis(2)


plot(pollock$Length,pollock$TotWt*24*0.0134*exp(.115*pollock$GearTemp)/pollock$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[1],ylim=c(0,.03),line=-1,,axes=FALSE);axis(1);axis(2)
points(S.1$Length,S.1$C,type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1)
lines(S.1$Length,S.1$C.hat,pch=16,col=col2,cex=1,lty=1,lwd=2)
S.a<-S.1[S.1$Length>5,]
S.a<-S.a[S.a$Length<78,]
model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT),logW=log(S.a$Wt)))
CA.a<-exp(coef(model.a)[1])
CB.a<-coef(model.a)[2]
lines(1:40,rep(.41/100,length(1:40)),col=col3,lwd=2,lty=1)
lines(40:130,rep(.46/100,length(40:130)),col=col3,lwd=2,lty=1)
#lines(S.1$Length,S.1$C,pch=16,col=colors()[300],cex=1,lty=1,lwd=2)
legend(45,.03,c("mean data (Stom.*digest)","VonB","bioenergetics","livingston","MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
lines(S.1$Length,S.1$C.MSM,col=col4,lwd=2)
#points(S.1$Length,S.1$C,type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1)

#Kerim's code
#Max Likelihood fit for cod
vonb<-list(Winf=21448.1, K=0.57325, t0=-0.495495, dexp=0.742355, H=7.48695)
# Important note:  "K" here is not the same value as the traditional
# vonB K; it's about 3 * higher (see Essington paper)
# H as consumption estimate is derived quantity:
# H = vonb$K * vonb$Winf^(1.-vonb$dexp)
#Plotting weight versus age
#  xage <- seq(0,20,0.1)
#  yage <- log (vonb$Winf) + (1./(1.- vonb$dexp)) * log ((1.-exp(-vonb$K*(1.-vonb$dexp)*(xage-vonb$t0))))
#plot(xage,yage)
# Indigestable stuff is 40%
 AAcons <- 0.6
# length/weight (cm to g) fit with same data
 lw_a <- 0.003932293 # mine =0.004117811
 lw_b <- 3.257113 #mine = 3.253258
lenlist      <- 1:130
wlist        <- lw_a*lenlist^lw_b
dailycons    <- (vonb$H * (wlist^vonb$dexp) / AAcons) / 365
dailyconsbio <- dailycons/wlist
#plot(lenlist,dailyconsbio)

plot(pcod$Length,pcod$TotWt*24*0.0134*exp(.115*pcod$GearTemp)/pcod$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",line=0,main="Pacific cod",xlim=c(0,130),ylim=c(0,.03),axes=FALSE);axis(1);axis(2)
points(S.2$Length[S.2$C<0.03],S.2$C[S.2$C<0.03],type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1)
lines(lenlist,dailyconsbio,pch=16,col=col2,cex=1,lty=1,lwd=2)
S.a<-S.2
S.a<-S.a[S.a$Length<98,]
model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT),logW=log(S.a$Wt)))
CA.a<-exp(coef(model.a)[1])
CB.a<-coef(model.a)[2]
lines(lenlist,0.03115127*(wlist^-0.1342909),lwd=2,col=col2,lty=2)  # data from Paul 
Pvals<-dailyconsbio/(0.03115127*(wlist^-0.1342909))
lines(1:55,rep(.62/100,length(1:55)),col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130)),col=col3,lwd=2,lty=1)
lines(S.2$Length,S.2$C.MSM,col=col4,lwd=2)


plot(arrowtooth$Length,arrowtooth$TotWt*24*0.0134*exp(.115*arrowtooth$GearTemp)/arrowtooth$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",line=-1,main=sp.dat$spp[3],ylim=c(0,.03),axes=FALSE);axis(1);axis(2)
points(S.3$Length,S.3$C,type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1)
S.a<-S.3[S.3$Length>5,]
S.a<-S.a[S.a$Length<78,]
model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT),logW=log(S.a$Wt)))
CA.a<-exp(coef(model.a)[1])
CB.a<-coef(model.a)[2]
lines(S.a$Length,sp.dat$CA[3]*(S.a$Wt^sp.dat$CB[3])*exp(3^0.3763),lwd=2,col=col2,lty=1)
lines(S.a$Length,sp.dat$CA[3]*(S.a$Wt^sp.dat$CB[3])*exp(6^0.3763),lwd=2,col=col2,lty=1)

lines(1:40,rep(.61/100,length(1:40)),col=col3,lwd=2,lty=1)
lines(40:130,rep(.63/100,length(40:130)),col=col3,lwd=2,lty=1)
#lines(S.3$Length,S.3$C,pch=16,col=colors()[300],cex=1,lty=1,lwd=2)
S.a<-S.3
S.a<-S.a[S.a$Length<78,]
lines(S.3$Length,S.3$C.MSM,col=col4,lwd=2)
#points(S.3$Length,S.3$C,type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1)

mtext(bquote(bold("Length (cm)")),side=1,line=1,outer=TRUE)
mtext(bquote(bold("Daily ration (g/g)")),side=2,line=1,outer=TRUE)
dev.print(device=postscript, paste("fig2ration.eps",sep=""), onefile=FALSE, horizontal=FALSE)

#####################################
## Fig 3 for Kerim - C in g/g and C in g for P cod only

# FIG 2 : plot of Ration and models (where data = mean S * digestion)
#graphics.off()

dev.new(width=5,height=6.5)
	layout(rbind(c(1,0),
				c(2,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(2)
col1<-colors()[351]  # data point colors
col2<-colors()[136] # red line for Cmax
col3<-colors()[496] # livingston lines
col4<-colors()[638]
#Kerim's code
#Max Likelihood fit for cod
vonb<-list(Winf=21448.1, K=0.57325, t0=-0.495495, dexp=0.742355, H=7.48695)
# Important note:  "K" here is not the same value as the traditional
# vonB K; it's about 3 * higher (see Essington paper)
# H as consumption estimate is derived quantity:
# H = vonb$K * vonb$Winf^(1.-vonb$dexp)
#Plotting weight versus age
#  xage <- seq(0,20,0.1)
#  yage <- log (vonb$Winf) + (1./(1.- vonb$dexp)) * log ((1.-exp(-vonb$K*(1.-vonb$dexp)*(xage-vonb$t0))))
#plot(xage,yage)
# Indigestable stuff is 40%
 AAcons <- 0.6
# length/weight (cm to g) fit with same data
 lw_a <- 0.003932293 # mine =0.004117811
 lw_b <- 3.257113 #mine = 3.253258
lenlist      <- 1:130
wlist        <- lw_a*lenlist^lw_b
dailycons    <- (vonb$H * (wlist^vonb$dexp) / AAcons) / 365
dailyconsbio <- dailycons/wlist
#plot(lenlist,dailyconsbio)

plot(pcod$Length,pcod$TotWt*24*0.0134*exp(.115*pcod$GearTemp)/pcod$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",line=0,main="Pacific cod",xlim=c(0,130),ylim=c(0,.03),axes=FALSE);axis(1);axis(2)
points(S.2$Length[S.2$C<0.03],S.2$C[S.2$C<0.03],type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1)
lines(lenlist,dailyconsbio,pch=16,col=col2,cex=1,lty=1,lwd=2)
#lines(S.2$Length,S.2$C.hat,pch=16,col=col2,cex=1,lty=1,lwd=2)
S.a<-S.2
S.a<-S.a[S.a$Length<98,]
model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT),logW=log(S.a$Wt)))
CA.a<-exp(coef(model.a)[1])
CB.a<-coef(model.a)[2]
#lines(S.a$Length,CA.a*(S.a$Wt^CB.a)*mean(Temp.data$fT1),lwd=2,col=col2,lty=2)
lines(lenlist,0.03115127*(wlist^-0.1342909),lwd=2,col=col2,lty=2)  # data from Paul 
Pvals<-dailyconsbio/(0.03115127*(wlist^-0.1342909))
lines(1:55,rep(.62/100,length(1:55)),col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130)),col=col3,lwd=2,lty=1)
lines(S.2$Length,S.2$C.MSM,col=col4,lwd=2)
#legend(45,.03,c("Mean observed","Cmax (from Lab)","C from Von B","Livingston","MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
mtext(bquote(bold("Daily ration (g/g)")),side=2,line=2,outer=FALSE)

plot(pcod$Length,pcod$TotWt*24*0.0134*exp(.115*pcod$GearTemp),pch=16,col=col1,cex=.5,ylab="",xlab="",line=1,main="",xlim=c(0,130),ylim=c(0,400),axes=FALSE);axis(1);axis(2)
points(S.2$Length[S.2$C<0.03],S.2$C[S.2$C<0.03]*S.2$Wt[S.2$C<0.03],type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1)
lines(lenlist,dailycons,pch=16,col=col2,cex=1,lty=1,lwd=2)
#lines(S.2$Length,S.2$C.hat,pch=16,col=col2,cex=1,lty=1,lwd=2)
S.a<-S.2
S.a<-S.a[S.a$Length<98,]
model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT),logW=log(S.a$Wt)))
CA.a<-exp(coef(model.a)[1])
CB.a<-coef(model.a)[2]
#lines(S.a$Length,CA.a*(S.a$Wt^CB.a)*mean(Temp.data$fT1),lwd=2,col=col2,lty=2)
lines(lenlist,0.03115127*(wlist^-0.1342909)*wlist,lwd=2,col=col2,lty=2)  # data from Paul 

Pvals<-dailyconsbio/(0.03115127*(wlist^-0.1342909))
lines(1:55,rep(.62/100,length(1:55))*wlist[1:55],col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130))*wlist[55:130],col=col3,lwd=2,lty=1)
lines(S.2$Length,S.2$C.MSM*S.2$Wt,col=col4,lwd=2)
legend(5,375,c("Mean observed","Cmax (from Lab)","C from Von B","Livingston","MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
mtext(bquote(bold("Daily ration (g)")),side=2,line=2,outer=FALSE)

mtext(bquote(bold("Length (cm)")),side=1,line=1,outer=TRUE)

dev.print(device=postscript, paste("fig4Kerim.eps",sep=""), onefile=FALSE, horizontal=FALSE)


##V2 - without points
#plot(S.2$Length[S.2$C<0.03],S.2$C[S.2$C<0.03],type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1,ylab="",xlab="",line=0,main="Pacific cod",xlim=c(0,130),ylim=c(0,.03),axes=FALSE);axis(1);axis(2)
#lines(lenlist,dailyconsbio,pch=16,col=col2,cex=1,lty=1,lwd=2)
##lines(S.2$Length,S.2$C.hat,pch=16,col=col2,cex=1,lty=1,lwd=2)
#S.a<-S.2
#S.a<-S.a[S.a$Length<98,]
#model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT),logW=log(S.a$Wt)))
#CA.a<-exp(coef(model.a)[1])
#CB.a<-coef(model.a)[2]
##lines(S.a$Length,CA.a*(S.a$Wt^CB.a)*mean(Temp.data$fT1),lwd=2,col=col2,lty=2)
#lines(lenlist,0.03115127*(wlist^-0.1342909),lwd=2,col=col2,lty=2)  # data from Paul 
#Pvals<-dailyconsbio/(0.03115127*(wlist^-0.1342909))
#lines(1:55,rep(.62/100,length(1:55)),col=col3,lwd=2,lty=1)
#lines(55:130,rep(1.33/100,length(55:130)),col=col3,lwd=2,lty=1)
#lines(S.2$Length,S.2$C.MSM,col=col4,lwd=2)
##legend(45,.03,c("Mean observed","Cmax (from Lab)","C from Von B","Livingston","MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
#mtext(bquote(bold("Daily ration (g/g)")),side=2,line=2,outer=FALSE)

#plot(S.2$Length[S.2$C<0.03],S.2$C[S.2$C<0.03]*S.2$Wt[S.2$C<0.03],type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1,ylab="",xlab="",line=1,main="",xlim=c(0,130),ylim=c(0,400),axes=FALSE);axis(1);axis(2)
#lines(lenlist,dailycons,pch=16,col=col2,cex=1,lty=1,lwd=2)
##lines(S.2$Length,S.2$C.hat,pch=16,col=col2,cex=1,lty=1,lwd=2)
#S.a<-S.2
#S.a<-S.a[S.a$Length<98,]
#model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT)),logW=log(S.a$Wt)))
#CA.a<-exp(coef(model.a)[1])
#CB.a<-coef(model.a)[2]
##lines(S.a$Length,CA.a*(S.a$Wt^CB.a)*mean(Temp.data$fT1),lwd=2,col=col2,lty=2)
#lines(lenlist,0.03115127*(wlist^-0.1342909)*wlist,lwd=2,col=col2,lty=2)  # data from Paul 
#Pvals<-dailyconsbio/(0.03115127*(wlist^-0.1342909))
#lines(1:55,rep(.62/100,length(1:55))*wlist[1:55],col=col3,lwd=2,lty=1)
#lines(55:130,rep(1.33/100,length(55:130))*wlist[55:130],col=col3,lwd=2,lty=1)
#lines(S.2$Length,S.2$C.MSM*S.2$Wt,col=col4,lwd=2)
#legend(5,375,c("Mean observed","Cmax (from Lab)","C from Von B","Livingston","MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
#mtext(bquote(bold("Daily ration (g)")),side=2,line=2,outer=FALSE)
#V2 - without points

plot(S.2$Length[S.2$C<0.03],S.2$C[S.2$C<0.03]*91.25,type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1,ylab="",xlab="",line=0,main="Pacific cod",xlim=c(0,130),ylim=c(0,3),axes=FALSE);axis(1);axis(2)
lines(lenlist,dailyconsbio*365,pch=16,col=col2,cex=1,lty=1,lwd=2)
#lines(S.2$Length,S.2$C.hat,pch=16,col=col2,cex=1,lty=1,lwd=2)
S.a<-S.2
S.a<-S.a[S.a$Length<98,]
model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT),logW=log(S.a$Wt)))
CA.a<-exp(coef(model.a)[1])
CB.a<-coef(model.a)[2]
#lines(S.a$Length,CA.a*(S.a$Wt^CB.a)*mean(Temp.data$fT1),lwd=2,col=col2,lty=2)
lines(lenlist,0.03115127*(wlist^-0.1342909)*91.25,lwd=2,col=col2,lty=2)  # data from Paul 
Pvals<-dailyconsbio/(0.03115127*(wlist^-0.1342909))
lines(1:55,rep(.62/100,length(1:55))*91.25,col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130))*91.25,col=col3,lwd=2,lty=1)
lines(S.2$Length,S.2$C.MSM*91.25,col=col4,lwd=2)
#legend(45,.03,c("Mean observed","Cmax (from Lab)","C from Von B","Livingston","MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
mtext(bquote(bold("Annual ration (g/g)")),side=2,line=2,outer=FALSE)


plot(S.2$Length[S.2$C<0.03],S.2$C[S.2$C<0.03]*S.2$Wt[S.2$C<0.03]*91.25/1000,type="b",pch=16,col=colors()[300],cex=.7,lty=1,lwd=1,ylab="",xlab="",line=1,main="",xlim=c(0,130),ylim=c(0,40),axes=FALSE);axis(1);axis(2)
lines(lenlist,dailycons*365/1000,pch=16,col=col2,cex=1,lty=1,lwd=2)
#lines(S.2$Length,S.2$C.hat,pch=16,col=col2,cex=1,lty=1,lwd=2)
S.a<-S.2
S.a<-S.a[S.a$Length<98,]
model.a<-lm(logC~logW,data=list(logC=log(S.a$C/S.a$fT),logW=log(S.a$Wt)))
CA.a<-exp(coef(model.a)[1])
CB.a<-coef(model.a)[2]
#lines(S.a$Length,CA.a*(S.a$Wt^CB.a)*mean(Temp.data$fT1),lwd=2,col=col2,lty=2)
lines(lenlist,0.03115127*(wlist^-0.1342909)*wlist*91.25/1000,lwd=2,col=col2,lty=2)  # data from Paul 
Pvals<-dailyconsbio/(0.03115127*(wlist^-0.1342909))
lines(1:55,rep(.62/100,length(1:55))*wlist[1:55]*91.25/1000,col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130))*wlist[55:130]*91.25/1000,col=col3,lwd=2,lty=1)
lines(S.2$Length,S.2$C.MSM*S.2$Wt*91.25/1000,col=col4,lwd=2)
legend(5,30,c("Mean observed","Cmax (from Lab)","C from Von B","Livingston","MSM"),col=c(colors()[300],col2,col2,col3,col4),lty=c(1,2,1,1,1),lwd=2,box.lwd=0,box.lty=0)
mtext(bquote(bold("Annual ration (Kg)")),side=2,line=2,outer=FALSE)
mtext(bquote(bold("Length (cm)")),side=1,line=1,outer=TRUE)

dev.print(device=postscript, paste("fig4Kerimv2.eps",sep=""), onefile=FALSE, horizontal=FALSE)




# FIG 3

dev.new(width=4,height=6.5)
par(mfrow=c(3,1))
ST.1<-ST.1[ST.1$T>-5,]
plot(ST.1$T,ST.1$S*0.0143*exp(0.115*ST.1$T)*24,pch=16,col="black",cex=1,ylab="",xlab="Gear Temp",main=sp.dat$spp[1],xlim=c(-3,11),axes=FALSE);axis(1);axis(2)
par(new=TRUE)
plot(-3:20,bioenergetics(par=model.parms[1,],W=S.1$Wt[1],TempC=-3:20,P=1,Eprey=3000,Efish=3000)$fTc,type="l",ylim=c(0.4,1.3),ylab="",xlab="",lwd=2,axes=FALSE,xlim=c(-3,11));axis(1);axis(4)

plot(ST.2$T,ST.2$S*0.0143*exp(0.115*ST.2$T)*24,pch=16,col="black",cex=1,ylab="max daily C (g/g)",xlab="Gear Temp",main=sp.dat$spp[2],xlim=c(-3,11),axes=FALSE);axis(1);axis(2)
par(new=TRUE)
plot(-3:20,bioenergetics(par=model.parms[4,],W=S.1$Wt[1],TempC=-3:20,P=1,Eprey=3000,Efish=3000)$fTc,type="l",ylim=c(0.4,1),ylab="",xlab="",lwd=2,axes=FALSE,xlim=c(-3,11));axis(1);axis(4)

plot(ST.3$T,ST.3$S*0.0143*exp(0.115*ST.3$T)*24,pch=16,col="black",cex=1,ylab="",xlab="Gear Temp",main=sp.dat$spp[3],xlim=c(-3,11),axes=FALSE);axis(1);axis(2)
par(new=TRUE)
plot(-3:20,bioenergetics(par=model.parms[5,],W=S.1$Wt[1],TempC=-3:20,P=1,Eprey=3000,Efish=3000)$fTc,type="l",ylim=c(0,1),ylab="",xlab="",lwd=2,axes=FALSE,xlim=c(-3,11));axis(1);axis(4)
dev.print(device=postscript, paste("fig3fT.eps",sep=""), onefile=FALSE, horizontal=FALSE)



#FIG. 4 Upaij as a function of size
dev.new(width=6,height=6.5)
	layout(rbind(c(1,0,2),
				c(3,0,4),
				c(5,0,6))
		,widths=c(1,0.01,1),
		heights=rep(1,3))
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(5,5,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(6)
	
	#dat1<-na.omit(data.frame(Length=as.numeric(names(test.11$wtot.ratio.mean)),plk.prop=test.11$wtot.ratio.mean,pcod.prop=test.12$wtot.ratio.mean))
	#m<-glm(cbind(plk.prop,pcod.prop,1-(plk.prop+pcod.prop))~length,data=list(length=dat1$Length,plk.prop=dat1$plk.prop,pcod.prop=dat1$pcod.prop),family=Gamma(link="log"))

for (predd in 1:3){
	for (preyy in 1:2){
		eval(parse(text=paste("dat.use<-test.",predd,preyy,"$wtot.ratio.mean",sep="")))
		dat1<-na.omit(data.frame(Length=as.numeric(names(dat.use)),prop=dat.use))
		#dat1<-dat1[dat1$prop>0,]
		dat1$prop[dat1$prop==0]<-0.0001
		if(preyy>1){
			dat1<-dat1[dat1$prop<0.04,]			
		}
		m3<-glm(prop~length,data=list(length=dat1$Length,prop=dat1$prop),family=Gamma(link="log"))
		#m3<-glm(cbind(prop,1-prop)~length,data=list(length=dat1$Length,prop=dat1$prop),family=binomial(link="log"))
		new<-data.frame(loglength=log(1:130),length=1:130,length2=(1:130)^2)
		new$prop.hat<-exp(predict(m3,newdata=new))
		new$prop.hat[new$prop.hat>max.na(dat.use)]<-max.na(dat.use)
		if(preyy>1){
			new$prop.hat[new$length>70]<-0	
		}
		if(predd==1){
			new$prop.hat[new$length>81]<-0	
		}
		plot(na.omit(data.frame(Length=as.numeric(names(dat.use)),prop=dat.use)),
				pch=16, line=-1,main=paste(sp.dat$spp[predd],sp.dat$spp[preyy],sep="-->"),ylab="",xlab="",axes=FALSE);axis(1);axis(2)
		points(new$length,new$prop.hat,col="blue",type="l",lwd=2)
		points(dat1$Length,fitted(m3),col="red",type="l",lwd=2)
	}
}
mtext(bquote(bold("Length (cm)")),side=1,line=2,outer=TRUE)
mtext(bquote(bold("Stomach proportion (by Wt.)")),side=2,line=2,outer=TRUE)

dev.print(device=postscript, paste("fig4K.eps",sep=""), onefile=FALSE, horizontal=FALSE)



dev.new(width=6,height=6.5)
	layout(rbind(c(1,0,2),
				c(3,0,4),
				c(5,0,6))
		,widths=c(1,0.01,1),
		heights=rep(1,3))
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(5,5,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(6)

for (predd in 1:3){
	for (preyy in 1:2){
		eval(parse(text=paste("dat.use<-test.",predd,preyy,"$wtot.ratio.mean",sep="")))
		dat1<-na.omit(data.frame(Length=as.numeric(names(dat.use)),prop=dat.use))
		dat1$prop[dat1$prop==0]<-0.0001
		if(preyy>1){
			dat1<-dat1[dat1$prop<0.04,]			
		}
		m3<-glm(prop~length,data=list(length=dat1$Length,prop=dat1$prop),family=Gamma(link="log"))
		new<-data.frame(loglength=log(1:130),length=1:130,length2=(1:130)^2)
		new$prop.hat<-exp(predict(m3,newdata=new))
		new$prop.hat[new$prop.hat>max.na(dat.use)]<-max.na(dat.use)
		if(preyy>1){
			new$prop.hat[new$length>70]<-0	
		}
		if(predd==1){
			new$prop.hat[new$length>81]<-0	
		}
		plot(density(dat1$prop), line=-1,main=paste(sp.dat$spp[predd],sp.dat$spp[preyy],sep="-->"),ylab="",xlab="",axes=FALSE);axis(1);axis(2)
		shape<-as.numeric(gamma.shape(m3)[1])
		scale<-as.numeric(mean(dat1$prop)/shape)
		rate<-1/scale
		disp<-1/shape
		points(density(rgamma(10000,1/disp,scale=mean(dat1$prop)*disp)),type="l",col="red")
		
		
	}
}
mtext(bquote(bold("Length (cm)")),side=1,line=2,outer=TRUE)
mtext(bquote(bold("Stomach proportion (by Wt.)")),side=2,line=2,outer=TRUE)
dev.print(device=postscript, paste("fig5Kdensity.eps",sep=""), onefile=FALSE, horizontal=FALSE)

## plot predicted Upaij.hat


##########################################
## Essingtion Von B ration using Length at age
##########################################
LA_data2<-data.frame(Ages=plk_WA$AGE,Lengths=plk_WA$LENGTH/10)
LA_data<-LA_data2[LA_data2$Ages<24,]
WA_data2<-data.frame(Ages=plk_WA$AGE,Weights=plk_WA$WEIGHT)
WA_data<-WA_data2[WA_data2$Ages<24,]
WA_data<-na.omit(WA_data)

rm(LA_data2)

Ustart<-tapply(LA_data$Lengths,LA_data$Ages,median)*beta
Ustart<-c(Ustart[-(length(Ustart))],10^10)
U<-Ustart
beta<-Ustart
X<-LA_data$Lengths
y<-LA_data$Ages
llik.oprobit <- function(par, X, y) {
	beta <- par
	Y <- as.matrix(y)
	X1 <- as.matrix(X)
	beta11 <- beta[1:11]
	score <- as.vector(X1%*%beta11)
	p1a <- pnorm(beta[12] - score)
	p2a <- pnorm(beta[13] - score)
	p3a <- pnorm(beta[14] - score)
	p4a <- pnorm(beta[15] - score)
	p1 <- log(p1a)
	p2 <- log(p2a - p1a)
	p3 <- log(p3a - p2a)
	p4 <- log(p4a - p3a)
	p5 <- log(1 - p4a)
	phi <- (Y==1)*p1 + (Y==2)*p2 + (Y==3)*p3 + (Y==4)*p4 + (Y==5)*p5
	return(-sum(phi))
}
#start <- c(lm(LA_data$Ages ~ LA_data$Lengths)$coefficients[-1],2,4,6,7)
#result <- optim(start, llik.oprobit, y=LA_data$Ages, X=LA_data$Lengths, method="BFGS", hessian=TRUE)

vonBL<-function(data,pars){
	# LAdata is a dataframe with Lengths and Ages are vectors of weight at age data
	# use ordered probit to fit regression
	LA_data<-data$LA_data
	age<-as.numeric(LA_data$Ages)
	Length.obs<-LA_data$Lengths
	a<-data$a
	b<-data$b
	H<-pars[1]
	K<-pars[2]
	t0<-pars[3]
	logsigma<-pars[4]
	sigma<-exp(logsigma)
	E<-(H/b)*a^(1/3)
	Linf<-E/K
	nobs<-length(LA_data$Lengths)
	Length.hat<-rep(0,nobs)
	LL<-rep(0,nobs)
	for (i in 1:nobs){
		Length.hat[i]<-Linf*(1-exp(-K*(age[i]-t0)))
		LL[i]<-dnorm(Length.hat[i]-Length.obs[i],sigma,log=TRUE)
	}
	-sum(LL)
}

vonB(data=list(LA_data=LA_data,a=sp.dat$LW.a[1],b=sp.dat$LW.b[1]),pars=c(1,1,0,log(.2)))
result<-optim(c(1,1,0,log(.2)),vonB,data=list(LA_data=LA_data,a=sp.dat$LW.a[1],b=sp.dat$LW.b[1]), method="BFGS", hessian=TRUE)

vonB.predict<-function(data,pars){
	# LAdata is a dataframe with Lengths and Ages are vectors of weight at age data
	# use ordered probit to fit regression
	LA_data<-data$LA_data
	A<-data$A
	d<-data$d
	age<-as.numeric(LA_data$Ages)
	Length.obs<-LA_data$Lengths
	a<-data$a
	b<-data$b
	H<-pars[1]
	k<-pars[2]
	K<-k/d
	t0<-pars[3]
	logsigma<-pars[4]
	sigma<-exp(logsigma)
	E<-(H/b)*a^(d-1)
	m<-d*b+(1-b)
	Linf<-(E/K)^(1/(1-m))
	nobs<-length(LA_data$Lengths)
	Length.hat<-rep(0,nobs)
	LL<-rep(0,nobs)
	for (i in 1:nobs){
		Length.hat[i]<-Linf*((1-exp(-K*(1-m)*(age[i]-t0)))^(1/(1-m)))
		LL[i]<-dnorm(Length.hat[i]-Length.obs[i],sigma,log=TRUE)
	}
	NLL<--sum(LL)
	W<-a*Length.obs^b
	W.hat<-a*Length.hat^b
	annualC<-(H/A)*W^d
	C.daily<-annualC/365
	return(list(Length.hat=Length.hat,NLL=NLL,Linf=Linf,W=W,W.hat=W.hat,C.daily=C.daily,annualC=annualC))
}


L.hat<-vonB.predict(data=list(LA_data=LA_data,a=sp.dat$LW.a[1],b=sp.dat$LW.b[1],A=0.65,d=2/3),pars=result$par)
plot(LA_data)
points(LA_data$Ages,L.hat$Length.hat,col="red",lwd=2,pch=16)
vonB.plk<-result$par
names(vonB.plk)<-c("H","k","logsigma")
vonB.plk<-data.frame(t(vonB.plk))
vonB.plk$Linf<-(vonB.plk$H/(sp.dat$LW.b[1]))*(sp.dat$LW.a[1])^(1/3)/(vonB.plk$k)

LA_data2<-data.frame(Ages=pcod_WA$AGE,Lengths=pcod_WA$LENGTH/10)
LA_data<-LA_data2[LA_data2$Ages<24,]
rm(LA_data2)
result<-optim(c(1,1,0,log(.2)),vonB,data=list(LA_data=LA_data,a=sp.dat$LW.a[2],b=sp.dat$LW.b[2]), method="BFGS", hessian=TRUE)
vonB.pcod<-result$par
names(vonB.pcod)<-c("H","k","logsigma")
vonB.pcod<-data.frame(t(vonB.pcod))
vonB.pcod$Linf<-(vonB.pcod$H/(sp.dat$LW.b[2]))*(sp.dat$LW.a[2])^(1/3)/(vonB.pcod$k)

#Kerim's code
#Max Likelihood fit for cod
vonb<-list(Winf=21448.1, K=0.57325, t0=-0.495495, dexp=0.742355, H=7.48695)


###########################################################################
## Final graph of rations
###########################################################################

dev.new(width=4,height=6.5)
	layout(rbind(c(1,0),
				c(2,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(2)
col1<-colors()[351]
col2<-colors()[136]
col3<-colors()[496]
col4<-colors()[638]
col5<-colors()[499]

## g/g first
## POLLOCK
WAdat<-plk_WA
LAdat<-plk_LA
LA_data2<-data.frame(Ages=dat$AGE,Lengths=dat$LENGTH/10)
LA_data<-LA_data2[LA_data2$Ages<24,]
WA_data2<-data.frame(Ages=dat$AGE,Weights=dat$WEIGHT)
WA_data<-WA_data2[WA_data2$Ages<24,]
WA_data<-na.omit(WA_data)
dat<-pollock
vonb<-reptoRlist("/Users/kkari/Documents/science/Projects/MSM/msm_data/vonb/s_vonbplk.rep")
Lengths<-0:100
a<-sp.dat$LW.a[1]; b<-sp.dat$LW.b[1]
Weights<-a*Lengths^b
H<-(vonb$H) #230.8486792
A<-0.65
d<-2/3
annualC<-((H/A)*(Weights)^d)/Weights
C<-annualC/365
C.obs<-dat$TotWt*24*0.0134*exp(.115*dat$GearTemp)/dat$Wt # g/g/d
C.obs.mean<-tapply(C.obs,dat$Length,mean.na)
Cmax.hat<-bioenergetics(model.parms[1,],a*(dat$Length^b),dat$GearTemp,1,5000,5000)$Cmax
Cmax.hat.mean<-tapply(Cmax.hat,dat$Length,mean.na)
Cmax.hat.juv<-bioenergetics(model.parms[2,],a*(dat$Length[dat$Length<30]^b),dat$GearTemp[dat$Length<30],1,5000,5000)$Cmax
Cmax.hat.juv.mean<-tapply(Cmax.hat.juv,dat$Length[dat$Length<30],mean.na)
Cmax.hat.adult<-bioenergetics(model.parms[2,],a*(dat$Length[dat$Length>=30]^b),dat$GearTemp[dat$Length>=30],1,5000,5000)$Cmax
Cmax.hat.adult.mean<-tapply(Cmax.hat.adult,dat$Length[dat$Length>=30],mean.na)
Y1<-bioenergetics(model.parms[1,],Weights,14,1,5000,5000)$Cmax/bioenergetics(model.parms[1,],Weights,14,1,5000,5000)$fTc
Y2.1<-bioenergetics(model.parms[2,],Weights[1:31],14,1,5000,5000)$Cmax/bioenergetics(model.parms[2,],Weights[1:31],14,1,5000,5000)$fTc
Y2.2<-bioenergetics(model.parms[2,],Weights[32:101],14,1,5000,5000)$Cmax/bioenergetics(model.parms[2,],Weights[32:101],14,1,5000,5000)$fTc
Y3<-C
plot(dat$Length,dat$TotWt*24*0.0134*exp(.115*dat$GearTemp)/dat$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[1],ylim=c(0,.03),line=-1,,axes=FALSE);axis(1);axis(2)
lines(as.numeric(names(C.obs.mean)),C.obs.mean)
lines(Lengths,Y1,col=col2,lwd=1,lty=2)
lines(as.numeric(names(Cmax.hat.mean)),Cmax.hat.mean,col=col2,lwd=2)
lines(as.numeric(names(Cmax.hat.juv.mean)),Cmax.hat.juv.mean,col=col4,lwd=2)
lines(as.numeric(names(Cmax.hat.adult.mean)),Cmax.hat.adult.mean,col=col4,lwd=2)
lines(Lengths[1:31],Y2.1,col=col4,lwd=1,lty=2)
lines(Lengths[32:101],Y2.2,col=col4,lwd=1,lty=2)
lines(as.numeric(names(Cmax.hat.mean)),Cmax.hat.mean*.3,col=col2,lwd=2,lty=3)
lines(Lengths,Y3,col=col5,lwd=2,lty=1)
lines(1:40,rep(.41/100,length(1:40)),col=col3,lwd=2,lty=1)
lines(40:130,rep(.46/100,length(40:130)),col=col3,lwd=2,lty=1)
legend(40,0.030,c("data","Ciannelli","Mason","vonB","livingston"),c("black",col2,col4,col4,col5,col3),lwd=2,lty=1,box.lwd=0)

## g/g first
# PCOD
WAdat<-pcod_WA
LAdat<-pcod_LA
mparms<-model.parms[2,]
LA_data2<-data.frame(Ages=dat$AGE,Lengths=dat$LENGTH/10)
LA_data<-LA_data2[LA_data2$Ages<24,]
WA_data2<-data.frame(Ages=dat$AGE,Weights=dat$WEIGHT)
WA_data<-WA_data2[WA_data2$Ages<24,]
WA_data<-na.omit(WA_data)
dat<-pcod
vonb<-reptoRlist("/Users/kkari/Documents/science/Projects/MSM/msm_data/vonb/s_vonbpcod.rep")
Lengths<-0:100
a<-sp.dat$LW.a[1]; b<-sp.dat$LW.b[1]
Weights<-a*Lengths^b
H<-log(230.8486792)
A<-0.65
d<-2/3
annualC<-((H/A)*(Weights)^d)/Weights
C<-annualC/365
C.obs<-dat$TotWt*24*0.0134*exp(.115*dat$GearTemp)/dat$Wt # g/g/d
C.obs.mean<-tapply(C.obs,dat$Length,mean.na)
Cmax.hat<-bioenergetics(model.parms[4,],a*(dat$Length^b),dat$GearTemp,1,5000,5000)$Cmax
Cmax.hat.mean<-tapply(Cmax.hat,dat$Length,mean.na)
Y1<-bioenergetics(model.parms[4,],Weights,14,1,5000,5000)$Cmax/bioenergetics(model.parms[1,],Weights,14,1,5000,5000)$fTc
Y3<-C
plot(dat$Length,dat$TotWt*24*0.0134*exp(.115*dat$GearTemp)/dat$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[2],ylim=c(0,.03),line=-1,,axes=FALSE);axis(1);axis(2)
lines(as.numeric(names(C.obs.mean)),C.obs.mean)
lines(Lengths,Y1,col=col2,lwd=1,lty=2)
lines(as.numeric(names(Cmax.hat.mean)),Cmax.hat.mean,col=col2,lwd=2)
lines(Lengths,Y3,col=col5,lwd=2,lty=1)
lines(1:55,rep(.62/100,length(1:55)),col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130)),col=col3,lwd=2,lty=1)

 AAcons <- 0.6
# length/weight (cm to g) fit with same data
 lw_a <- 0.003932293 # mine =0.004117811
 lw_b <- 3.257113 #mine = 3.253258
lenlist      <- 1:130
wlist        <- lw_a*lenlist^lw_b
dailycons    <- (vonb$H * (wlist^(d)) / AAcons) / 365
dailyconsbio <- dailycons/wlist
Y3<-dailyconsbio
lines(lenlist,Y3,col=col5,lwd=2,lty=1)

lines(1:55,rep(.62/100,length(1:55))*wlist[1:55],col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130))*wlist[55:130],col=col3,lwd=2,lty=1)



dev.new(width=4,height=6.5)
	layout(rbind(c(1,0),
				c(2,0))
		,widths=c(1,0.01),
		heights=1)
	par(mar=c(1,1,1,1)) # margins of graph (bottom,left, top, right)
	par(mgp=c(3,.5,0)) # axis margins, distance of lab, distance of numbers, axis line
	par(oma=c(3,3,1,0))# outer margins of graph (bottom,left, top, right)
	layout.show(2)
col1<-colors()[351]
col2<-colors()[136]
col3<-colors()[496]
col4<-colors()[638]
col5<-colors()[499]

## g first
## POLLOCK
WAdat<-plk_WA
LAdat<-plk_LA
LA_data2<-data.frame(Ages=dat$AGE,Lengths=dat$LENGTH/10)
LA_data<-LA_data2[LA_data2$Ages<24,]
WA_data2<-data.frame(Ages=dat$AGE,Weights=dat$WEIGHT)
WA_data<-WA_data2[WA_data2$Ages<24,]
WA_data<-na.omit(WA_data)
dat<-pollock
vonb<-reptoRlist("/Users/kkari/Documents/science/Projects/MSM/msm_data/vonb/s_vonbplk.rep")
Lengths<-0:100
a<-sp.dat$LW.a[1]; b<-sp.dat$LW.b[1]
Weights<-a*Lengths^b
H<-(vonb$H) #230.8486792
A<-0.65
d<-2/3
annualC<-((H/A)*(Weights)^d)
C<-annualC/365
C.obs<-dat$TotWt*24*0.0134*91.25*exp(.115*dat$GearTemp) # g
C.obs.mean<-tapply(C.obs,dat$Length,mean.na)
Cmax.hat<-(bioenergetics(model.parms[1,],a*(dat$Length^b),dat$GearTemp,1,5000,5000)$Cmax)*a*(dat$Length^b)
Cmax.hat.mean<-tapply(Cmax.hat,dat$Length,mean.na)
Cmax.hat.juv<-bioenergetics(model.parms[2,],a*(dat$Length[dat$Length<30]^b),dat$GearTemp[dat$Length<30],1,5000,5000)$Cmax*a*(dat$Length^b)
Cmax.hat.juv.mean<-tapply(Cmax.hat.juv,dat$Length[dat$Length<30],mean.na)
Cmax.hat.adult<-bioenergetics(model.parms[2,],a*(dat$Length[dat$Length>=30]^b),dat$GearTemp[dat$Length>=30],1,5000,5000)$Cmax*a*(dat$Length^b)
Cmax.hat.adult.mean<-tapply(Cmax.hat.adult,dat$Length[dat$Length>=30],mean.na)
Y1<-bioenergetics(model.parms[1,],Weights,14,1,5000,5000)$Cmax/bioenergetics(model.parms[1,],Weights,14,1,5000,5000)$fTc

Y2.1<-bioenergetics(model.parms[2,],Weights[1:31],14,1,5000,5000)$Cmax/bioenergetics(model.parms[2,],Weights[1:31],14,1,5000,5000)$fTc
Y2.2<-bioenergetics(model.parms[2,],Weights[32:101],14,1,5000,5000)$Cmax/bioenergetics(model.parms[2,],Weights[32:101],14,1,5000,5000)$fTc
Y3<-annualC
Y1<-Y1*Weights*91.25
Y2.1<-Y1*Weights[1:31]*91.25
Y2.2<-Y1*Weights[32:101]*91.25

plot(dat$Length,dat$TotWt*24*0.0134*91.25*exp(.115*dat$GearTemp),pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[1],ylim=c(0,20*90),line=-1,,axes=FALSE);axis(1);axis(2)
lines(as.numeric(names(C.obs.mean)),C.obs.mean)
lines(Lengths,Y1,col=col2,lwd=1,lty=2)
lines(as.numeric(names(Cmax.hat.mean)),Cmax.hat.mean,col=col2,lwd=2)
lines(as.numeric(names(Cmax.hat.juv.mean)),Cmax.hat.juv.mean,col=col4,lwd=2)
lines(as.numeric(names(Cmax.hat.adult.mean)),Cmax.hat.adult.mean,col=col4,lwd=2)
lines(Lengths[1:31],Y2.1*a*(dat$Length^b),col=col4,lwd=1,lty=2)
lines(Lengths[32:101],Y2.2*a*(dat$Length^b),col=col4,lwd=1,lty=2)
lines(as.numeric(names(Cmax.hat.mean)),Cmax.hat.mean*.3,col=col2,lwd=2,lty=3)
lines(Lengths,Y3,col=col5,lwd=2,lty=1)
lines(1:40,a*((1:40)^b)*(rep(.41/100,length(1:40))),col=col3,lwd=2,lty=1)
lines(40:130,rep(.46/100,length(40:130))*a*((40:130)^b),col=col3,lwd=2,lty=1)


## g/g first
# PCOD
WAdat<-pcod_WA
LAdat<-pcod_LA
mparms<-model.parms[2,]
LA_data2<-data.frame(Ages=dat$AGE,Lengths=dat$LENGTH/10)
LA_data<-LA_data2[LA_data2$Ages<24,]
WA_data2<-data.frame(Ages=dat$AGE,Weights=dat$WEIGHT)
WA_data<-WA_data2[WA_data2$Ages<24,]
WA_data<-na.omit(WA_data)
dat<-pcod
vonb<-reptoRlist("/Users/kkari/Documents/science/Projects/MSM/msm_data/vonb/s_vonbpcod.rep")
Lengths<-0:100
a<-sp.dat$LW.a[1]; b<-sp.dat$LW.b[1]
Weights<-a*Lengths^b
H<-log(230.8486792)
A<-0.65
d<-2/3
annualC<-((H/A)*(Weights)^d)/Weights
C<-annualC/365
C.obs<-dat$TotWt*24*0.0134*exp(.115*dat$GearTemp)/dat$Wt # g/g/d
C.obs.mean<-tapply(C.obs,dat$Length,mean.na)
Cmax.hat<-bioenergetics(model.parms[4,],a*(dat$Length^b),dat$GearTemp,1,5000,5000)$Cmax
Cmax.hat.mean<-tapply(Cmax.hat,dat$Length,mean.na)
Y1<-bioenergetics(model.parms[4,],Weights,14,1,5000,5000)$Cmax/bioenergetics(model.parms[1,],Weights,14,1,5000,5000)$fTc
Y3<-C
plot(dat$Length,dat$TotWt*24*0.0134*exp(.115*dat$GearTemp)/dat$Wt,pch=16,col=col1,cex=.5,ylab="",xlab="",main=sp.dat$spp[2],ylim=c(0,.03),line=-1,,axes=FALSE);axis(1);axis(2)
lines(as.numeric(names(C.obs.mean)),C.obs.mean)
lines(Lengths,Y1,col=col2,lwd=1,lty=2)
lines(as.numeric(names(Cmax.hat.mean)),Cmax.hat.mean,col=col2,lwd=2)
lines(Lengths,Y3,col=col5,lwd=2,lty=1)
lines(1:55,rep(.62/100,length(1:55)),col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130)),col=col3,lwd=2,lty=1)

 AAcons <- 0.6
# length/weight (cm to g) fit with same data
 lw_a <- 0.003932293 # mine =0.004117811
 lw_b <- 3.257113 #mine = 3.253258
lenlist      <- 1:130
wlist        <- lw_a*lenlist^lw_b
dailycons    <- (vonb$H * (wlist^(d)) / AAcons) / 365
dailyconsbio <- dailycons/wlist
Y3<-dailyconsbio
lines(lenlist,Y3,col=col5,lwd=2,lty=1)

lines(1:55,rep(.62/100,length(1:55))*wlist[1:55],col=col3,lwd=2,lty=1)
lines(55:130,rep(1.33/100,length(55:130))*wlist[55:130],col=col3,lwd=2,lty=1)





outfile<-"/Users/kkari/Documents/science/Projects/MSM/vonB/vonB.dat"
# K regression coef for each sppXspp combo where K <-a*pred_size^b
cat("# Nobs ",file=outfile,append=FALSE,sep="\n")
cat(length(LA_data[,1]),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Lengths ",file=outfile,append=TRUE,sep="\n")
cat((LA_data$Lengths),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Ages ",file=outfile,append=TRUE,sep="\n")
cat((LA_data$Ages),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")



WA_data2<-data.frame(Ages=plk_WA$AGE,Weights=plk_WA$WEIGHT)
WA_data<-WA_data2[WA_data2$Ages<24,]
WA_data<-na.omit(WA_data)
outfile<-"/Users/kkari/Documents/science/Projects/MSM/vonB/vonBwplk.dat"
# K regression coef for each sppXspp combo where K <-a*pred_size^b
cat("# Nobs ",file=outfile,append=FALSE,sep="\n")
cat(length(WA_data[,1]),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Weights ",file=outfile,append=TRUE,sep="\n")
cat((WA_data$Weights),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Ages ",file=outfile,append=TRUE,sep="\n")
cat((WA_data$Ages),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")


WA_data2<-data.frame(Ages=pcod_WA$AGE,Weights=pcod_WA$WEIGHT)
WA_data<-WA_data2[WA_data2$Ages<24,]
WA_data<-na.omit(WA_data)
outfile<-"/Users/kkari/Documents/science/Projects/MSM/vonB/vonBwpcod.dat"
# K regression coef for each sppXspp combo where K <-a*pred_size^b
cat("# Nobs ",file=outfile,append=FALSE,sep="\n")
cat(length(WA_data[,1]),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Weights ",file=outfile,append=TRUE,sep="\n")
cat((WA_data$Weights),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Ages ",file=outfile,append=TRUE,sep="\n")
cat((WA_data$Ages),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")

WA_data2<-data.frame(Ages=pcod_WA$AGE,Weights=pcod_WA$WEIGHT)
WA_data<-WA_data2[WA_data2$Ages<24,]
WA_data<-na.omit(WA_data)
outfile<-"/Users/kkari/Documents/science/Projects/MSM/vonB/vonBwatf.dat"
# K regression coef for each sppXspp combo where K <-a*pred_size^b
cat("# Nobs ",file=outfile,append=FALSE,sep="\n")
cat(length(WA_data[,1]),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Weights ",file=outfile,append=TRUE,sep="\n")
cat((WA_data$Weights),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")
cat("# Ages ",file=outfile,append=TRUE,sep="\n")
cat((WA_data$Ages),file=outfile,append=TRUE,sep=" ");cat("",file=outfile,append=TRUE,sep="\n")



